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Abstract 

The  use  of  control  variates  is  a  well-known  variance  reduction  technique  for  discrete 
event  simulation  experiments.  Currently,  internal  control  variates  are  used  almost 
exclusively  by  practitioners  and  researchers  when  using  control  variates.  The  primary 
objective  of  this  study  is  to  explore  the  variance  reduction  achieved  by  the  replicative  use 
of  an  external,  analytical  model  to  generate  control  variates.  Performance  for  the 
analytical  control  variates  is  compared  to  the  performance  of  typical  internal  and  external 
control  variates  for  both  an  open  and  a  closed  queueing  network.  Performance  measures 
used  are  confidence  interval  width  reduction,  realized  coverage,  and  estimated  Mean 
Square  Error.  Results  of  this  study  indicate  analytical  control  variates  achieve  comparable 
confidence  interval  width  reduction  with  internal  and  external  control  variates.  However, 
the  analytical  control  variates  exhibit  greater  levels  of  estimated  bias.  Possible  causes  and 
remedies  for  the  observed  bias  are  discussed  and  areas  for  future  research  and  use  of 
analytical  control  variates  conclude  the  study. 
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REPLICATIVE  USE  OF  AN  EXTERNAL  MODEL 
IN  SIMULATION  VARIANCE  REDUCTION 

1.  Introduction 


1. 1  Statement  of  the  Problem 

Air  Force  and  industry  analysts  often  use  computer  simulations  to  study  large  scale 
systems.  Decision  makers  use  the  results  of  these  simulations  to  set  policy  and  allocate 
scarce  resources.  In  many  applications,  analysts  focus  on  one  measurement  of 
effectiveness  as  the  basis  for  evaluating  competing  policies  or  allocations.  Obviously,  the 
better  an  analyst  can  estimate  this  measure,  the  better  the  final  decision. 

Simulations  produce  random  outputs,  not  exact  answers.  To  interpret  a  simulation 
output,  analysts  typically  use  statistical  confidence  intervals  based  on  the  mean  and 
variance  of  the  output  random  variable  in  order  to  estimate  its  true  expected  value. 
Analysts  may  not  have  the  time  or  budget  to  complete  enough  replications  to  produce  a 
sufficiently  small  confidence  interval.  By  using  an  appropriate  variance  reduction 
technique,  the  analyst  can  significantly  shrink  the  size  of  the  confidence  interval  or  reduce 
the  computational  effort  needed  to  obtain  an  estimate  of  specified  accuracy. 

Control  variates  comprise  one  method  that  has  been  shown  to  significantly  reduce 
variance.  Control  variates  are  usually  classified  as  either  internal  or  external  control 
variates.  Both  methods  attempt  to  take  advantage  of  correlation  between  selected  random 
variables  to  reduce  variance.  Internal  control  variates  can  be  the  input  random  variables, 
or  simple  functions  of  them.  Unfortunately,  selecting  a  subset  of  all  the  possible  internal 
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control  variates  for  a  given  simulation  study  that  produces  the  small  and  reliable 
confidence  intervals  can  be  a  difficult  problem  [14],  External  control  variates  require  the 
creation  of  a  simplified  analytical  model  of  the  system  under  study  and  an  associated 
simulation  model  based  on  the  analytical  model.  Applying  external  control  variates  to 
large  simulations  can  be  quite  difficult  since  common  random  numbers  must  be  applied  to 
each  simulation  [13],  For  these  reasons,  external  control  variates  are  rarely  used  in 
practical  applications  [17], 

By  using  an  external,  analytical  model  to  generate  a  single  control  variate  for  each 
replication  of  a  simulation,  it  may  be  possible  to  avoid  both  the  undesirable  impact  of 
using  multiple  control  variates  and  the  difficult  application  of  common  random  numbers. 
Sharon  [19]  first  used  Jackson  networks  to  generate  control  variates  for  queueing  network 
simulations.  Tomick,  Litko,  and  Bauer  [21]  demonstrated  that  in  small  queueing  networks 
significant  variance  reduction  can  be  achieved  by  generating  control  variates  using  an 
anal54ical  model.  Dietz  and  Harmonosky  [6]  achieved  significant  variance  reduction  in  an 
aircraft  sortie  generation  simulation  by  using  a  composite  control  variate  representing  a 
total  expected  “hands-on”  time  per  sortie.  If  this  same  technique  can  reduce  variance  in 
large  real  world  simulations,  study  times  can  be  significantly  reduced.  This  study  will 
explore  the  variance  reduction  achieved  with  the  replicative  use  of  an  external  model  to 
generate  control  variates  when  applied  to  simulations  of  varying  sizes.  The  variance 
reduction  achieved  using  this  method  will  be  compared  to  the  variance  reduction  obtained 
in  the  same  simulations  using  typical  internal  and  external  control  variate  techniques. 
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1.2  Background 

The  following  background  information  is  provided  as  motivation  for  the  development 
of  this  study.  First,  the  theory  of  control  variates  (particularly  how  they  reduce  variance 
and  how  simulation  confidence  intervals  are  formed)  is  presented,  followed  by  a  discussion 
of  internal  and  external  control  variates.  The  theoretical  basis  for  using  an  external  model 
to  generate  a  single  control  variate  is  discussed  next. 

1.2.1  Control  Variates.  This  variance  reduction  technique  attempts  to  take 
advantage  of  the  correlation  between  different  random  variables  in  order  to  obtain  a 
tighter  estimate  of  a  response  variable.  For  example,  suppose  an  analyst  is  trying  to 
estimate  ^  =  E[Y\,  where  the  random  variable  Y  is  the  waiting  time  for  a  customer  in  a 
bank  line.  To  reduce  the  variance,  the  analyst  can  then  use  another  random  variable,  say 
C,  that  the  analyst  believes  to  be  correlated  to  Y  and  that  has  a  known  expectation  (mz  = 
£[C].  A  new  random  variable  Y(b)  =  Y-  b(C  -  fic)  can  be  constructed  for  each  replication 
where 

£[r(i)]  =  £[r]-i(£[C]-^c)  (1.1) 

is  an  unbiased  estimator  for  //  for  any  real  number  b.  Then  since, 

Var{Y(b))  =  Var{Y)  +  b^Var{C)  -  2bCov{Y,  C)  ( 1 .2) 

Var{Y{b))  will  have  a  smaller  value  than  Far(Y)  if  and  only  if 

2aCov(Y,Q>b^Var(Q  (1.3) 
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which  will  be  the  case  if  7  and  C  are  highly  correlated  and  b  is  selected  appropriately.  Let 
P  be  the  optimal  value  for  h  that  minimizes  Var(Y(b)).  By  setting  the  derivative  of 
equation  (1.2),  (with  respect  to  b),  equal  to  zero  we  find 


Cov(Y,C) 

Var{Y) 


(1.4) 


Normally  Cov{Y,C)  and  Var(Y)  are  not  known.  Then  for  n  replications,  can  be 
estimated  by 


(1.5) 


y=i 

where  Y  is  the  mean  of  the  n  observations  of  Y  and  C  is  the  mean  of  the  n  observations 
of  C.  Then  for  the  same  n  replications  we  can  estimate  fx  with  the  controlled  estimate 


Y(h  =  Y-p(C-fXc)  (16) 

This  method  can  be  generalized  for  several  control  variates,  Ci,  C2,  ...,  Q  with 
respective  means  fiq  to 


W  =  l'-I*.(C.-/i.)  (1.7) 

k=\ 

where  the  ^i’s  are  any  real  numbers.  If  the  Z>*’s  are  estimated  in  the  same  manner  as 
equation  (1.6),  the  variance  of  Y0)  is  given  by  [14] 
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and  the  associated  confidence  interval  is 


W  ±  -  R‘W^[r]  (1.9) 

where  is  the  coefficient  of  correlation  between  Y  and  C  =  [Ci,  Q],  n  is  the  number 
of  replications,  and  q  is  the  number  of  control  variates.  (Note  equations  (1.8)  and  (1.9) 
are  not  computational  formulas.  The  computational  formulas  are  given  in  chapter  3.)  As 
each  additional  control  variate  is  added,  the  /-statistic  loses  an  additional  degree  of 
freedom  and  the  loss  term  (n-2)l(n-q-l)  increases,  while  may  increase.  Due  to  these 
countering  effects,  at  some  point  the  addition  of  another  control  variate  can  cause  the 
confidence  interval  to  grow  rather  than  decrease.  [13] 

1.2.2  Internal  and  External  Control  Variates.  Control  variates  can  be  obtained 
either  internally  or  externally.  Internal  control  variates  arise  from  the  simulation  itself  and 
are  normally  some  combination  of  the  input  random  variables  or  simple  functions  of  them. 
Examples  are  interarrival  times,  service  times,  and  demand  levels  since  their  expectations 
are  known  and  analysis  of  the  system  should  provide  the  sign  of  the  correlation  between 
them  and  the  output  random  variable.  As  demonstrated  in  the  previous  section,  as  more 
control  variates  are  added  to  induce  more  correlation,  the  corresponding  confidence 
interval  will  begin  to  grow  at  some  point.  Selection  of  the  best  subset  of  all  possible 
internal  control  variates  then  becomes  a  problem  of  its  own  [12]. 

External  control  variates  are  normally  obtained  by  simultaneously  simulating  a  similar 
system  that  has  an  analytical  solution.  By  using  common  random  numbers  the  output 
variable  of  the  simple  simulation  should  be  positively  correlated  to  the  output  variable  of 
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the  system  under  study  (hopefully  highly  correlated).  Using  external  control  variates  is 
usually  not  “cheap”  due  to  the  extra  work  required  to  develop  a  similar,  simpler  system, 
the  effort  required  to  synchronize  the  common  random  numbers,  and  the  additional 
computer  time  required  to  run  the  external  simulation  [6].  However,  it  can  pay  off  if  the 
resulting  variance  reduction  is  larger  than  that  of  the  internal  control  variate. 

1.2.3  External  Model  (Analytical)  Control  Variates.  Since  finding  the  “best”  subset 
of  internal  control  variates  or  synchronizing  common  random  can  be  difficult  tasks  in  large 
simulations,  control  variates  are  often  not  considered.  However,  Tomick  et  al.  [21]  and 
Dietz  and  Harmonosky  [6]  both  demonstrated  that  control  variates  could  be  obtained 
directly  from  an  analytical  model.  The  control  variates  obtained  in  each  paper  can  be 
viewed  simply  as  new  random  variables  that  are  functions  of  the  input  random  variables. 
This  type  of  control  variate  can  be  considered  a  hybrid  of  both  internal  and  external 
control  variates  since  it  is  found  using  an  external  model  --  like  external  control  variates  - 
and  is  a  function  of  the  input  random  variables  that  doesn’t  require  another  simulation 
model  —  like  internal  control  variates.  In  order  to  avoid  confusion,  control  variates  found 
in  this  manner  will  be  referred  to  as  analytical  control  variates  throughout  the  rest  of  this 
study. 

To  illustrate  how  an  analytical  control  variate  can  be  generated  consider  a  simulation 
model  where  Y  is  the  output  random  variable  of  interest  with  E{Y\  =  //.  Also,  let  X  =  (Xi, 
X2,  ...,  Xm)  be  the  vector  of  input  random  variables  that  drive  the  simulation  model  with 
E\X\  =  ftz , /^).  An  appropriate  analytical  model  is  essentially  a  function  of  some 
subset  of  the  same  input  random  variables  that  drive  the  simulation,  ideally  all  of  the  input 
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random  variables.  This  function  should  then  provide  an  approximation  of  the  output  of 
the  simulation.  To  illustrate,  represent  the  analytical  model  as  the  function  J{X)  =  Z.  Let 
==  =  Then  a  new  random  variable,  the  analytical  control  variate, 

can  be  constructed  as  before  by  Y(b)  =  Y-  b(J  (A!")  -  Z^) ,  so  that 

£[}'(*)]  =  £[K]-S(£[/(X)]-2,)  (l.IO) 

is  an  unbiased  estimator  of  /u  for  any  real  number  b. 

To  generate  analytical  control  variate  point  estimates  let  Xj  be  the  mean  of  the 
observed  values  of  X  for  the y-th  replication.  Further  let  f{Xj^  =  Zj  be  the  value  of  the 

analytical  model  for  the y-th  replication  and  let  Z  be  the  mean  of  these  calculations  for  all 
n  replications  performed.  Then  an  analytically  controlled  estimate  of  fx  is 

=  (1.11) 

where  p  is  estimated  using  equation  (1.5).  As  long  as  /(^/)  Xi  are  correlated 
variance  reduction  should  occur. 

1.3  Objective 

The  primary  objective  of  this  study  is  to  investigate  the  variance  reduction  of  the 
replicative  use  an  external  model  to  generate  analytical  control  variates  for  simulations. 
This  study  will  focus  on  simulations  with  a  single  random  variable  output.  The  primary 
performance  characteristic  for  this  study  will  be  confidence  interval  width  rather  than 
variance.  Confidence  interval  width  is  chosen  for  two  reasons.  First,  analysts  are 
primarily  concerned  with  confidence  interval  width  when  performing  simulation  studies. 
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Secondly,  confidence  interval  width  is  a  function  of  variance,  the  number  of  control 
variates  chosen,  and  number  of  replications.  Thus,  internal  control  variates  can  be 
compared  to  external  and  analytical  control  variates  in  a  meaningful  manner.  Comparisons 
of  the  point  estimate  confidence  interval  width  will  be  made  when  no  control  variates  are 
used  to  the  point  estimate  confidence  interval  width  when  internal  control  variates, 
external  control  variates,  and  analytical  control  variates  are  used  on  the  same  simulation. 
In  addition,  realized  confidence  interval  coverage  and  estimated  mean  square  error  will  be 
compared  for  all  cases  to  investigate  the  possibility  of  bias  error. 

1.3.1  Approach.  The  objectives  of  this  study  will  be  accomplished  in  the  following 
manner. 

1.  Create  a  simple  open  queueing  network  and  an  appropriate  simulation  model. 
Develop  an  external  analytical  model  whose  output  is  correlated  to  the  output  of 
the  simulation  model.  Create  a  separate  simulation  model  of  the  analytical  model. 

2.  Perform  simulation  experiments  on  three  different  design  points  of  the  open 
queuing  network  and  generate  internal,  analytical,  and  external  control  variates. 
Calculate  response  statistics  without  controls  and  with  controls  and  compare  the 
resulting  confidence  intervals. 

3.  Create  or  find  a  more  complex  closed  queueing  network  and  build  an  appropriate 
simulation  model.  Using  the  mean  value  analysis  algorithm  as  the  analytical  model, 
create  a  simulation  program  of  the  analytical  model. 

4.  Perform  simulation  experiments  on  twelve  different  design  points  of  the  closed 
queuing  network  and  generate  internal,  analytical,  and  external  control  variates. 
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Calculate  response  statistics  without  controls  and  with  controls  and  compare  the 
resulting  confidence  intervals. 

1.3.2  Scope.  Only  estimates  of  single  random  variable  confidence  intervals  will  be 
explored.  Estimating  the  reduction  in  multivariate  confidence  intervals  will  not  be 
discussed.  Estimating  simulation  outputs  can  be  accomplished  using  replication  methods, 
batch  means,  or  the  regenerative  method.  Only  the  replication  method  will  be  considered 
in  this  study.  Other  variance  reduction  techniques  such  as  common  random  numbers, 
antithetic  variates,  indirect  estimation,  or  conditional  expectations  will  not  be  considered. 
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2.  Previous  Work 


2.1  Introduction 

A  review  of  the  current  literature  on  the  use  of  control  variates  to  reduce  the  variance 
of  simulation  studies  as  it  pertains  to  this  study  is  presented  below.  The  discussion  will 
focus  on  the  use  of  control  variates  when  estimating  single  variate  responses  using 
replicative  simulation  methods.  The  first  section  presents  some  of  the  basic  references  on 
control  variates,  followed  by  a  review  on  the  use  of  external  control  variates.  Selected 
topics  on  the  use  of  specific  internal  control  variates  as  applicable  to  this  study  are 
presented  next.  The  final  section  reviews  the  literature  on  analytical  control  variates 

2.2  General  Control  Variate  Methodology 

2.2.1  Law  and  Kelton  (1991).  Several  authors  present  excellent  references  on  the 
use  of  control  variates  to  reduce  variance.  The  Law  and  Kelton  textbook  (1991)  on 
simulation  modeling  includes  a  section  on  the  use  of  control  variates  in  their  chapter  of 
variance  reduction  techniques.  They  provide  an  overview  on  control  variate  theory, 
application,  and  a  brief  discussion  on  the  difference  between  internal  and  external  control 
variates.  In  addition  to  their  excellent  coverage  of  control  variates,  they  also  provide  an 
extensive  reference  list  of  control  variate  literature. 

2.2.2  Nelson  (1987).  This  article  presents  a  guide  for  simulation  practitioners  for 
applying  three  variance  reduction  techniques,  including  control  variates.  Nelson  provides 
methods  for  finding  point  and  interval  estimators,  software  requirements,  and  guidelines 
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for  experiment  design  when  using  control  variates.  This  article  is  designed  as  a  tutorial, 
not  as  a  definitive  presentation  of  control  variate  theory.  As  such  it  is  very  useful  in 
understanding  how  the  theory  of  control  variates  is  applied  to  actual  simulation  studies. 
For  example,  Nelson  points  out  that  the  computation  of  control  variate  point  and  interval 
estimators  is  the  same  as  estimating  the  intercept  term  of  a  least-squares  regression  of  the 
responses  on  the  control  variates  minus  their  known  means. 

Nelson  goes  on  to  describe  two  methods  for  deciding  which  subset  of  possible 
control  variates  to  use  for  reducing  variance.  (It  should  be  noted  that  when  Nelson  refers 
to  control  variates,  he  is  referring  to  only  internal  control  variates.)  One  method  is  to  use 
a  regression  software  package  and  perform  stepwise  regression  on  all  possible  control 
variates  and  then  select  the  subset  of  controls  that  create  the  largest  amount  of  variance 
reduction.  He  also  proposes  a  means  of  determining  a  marginal  improvement  ratio  for 
adding  an  additional  control  variate  to  the  set  of  controls  already  in  use.  For  example,  for 
a  set  of  q  control  variates,  the  marginal  improvement  ratio  is  computed  by  comparing  1- 
for  the  set  of  ^  +  1  control  variates,  and  1  - ior  q  control  variates,  where  is  an 
estimate  of  the  square  of  the  multiple  correlation  coefficient  of  the  response  variable.  He 
then  provides  a  table  of  marginal  improvement  ratios  necessary  for  adding  an  additional 
control  variate  based  on  the  number  of  replications  performed  and  number  of  control 
variates  already  included. 

2.2.3  Lavenherg  and  Welch  (1981).  This  is  a  very  detailed  survey  on  the  application 
of  control  variates.  The  methods  they  present  in  generating  confidence  intervals  form  the 
foundation  of  the  methods  used  in  this  study.  As  in  Nelson  (1987),  the  focus  is  on  internal 
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control  variates,  but  Lavenberg  and  Welch  do  point  out  that  control  variates  can  be 
obtained  from  the  simulation  of  a  second  related  model  that  has  a  known  expected  value 
and  is  driven  by  the  same  sequences  of  random  numbers  that  are  used  in  the  simulation 
under  study.  These  control  variates  can  then  be  looked  upon  as  complex  functions  of  the 
input  random  variables  and  the  theory  of  control  variates  they  present  holds  for  these 
types  of  control  variates  too. 

Other  areas  investigated  include  techniques  for  generating  control  variates  (internal), 
inefficiencies  resulting  from  estimating  control  variate  coefficients,  a  review  of  the  studies 
completed  at  the  time  of  the  article,  and  directions  for  future  research.  Particularly  useful 
are  the  two  appendices  of  this  article  which  provide  detailed  discussions  on  the  application 
of  control  variates  to  generate  confidence  intervals.  The  first  appendix  describes  the 
generation  of  confidence  intervals  using  the  method  of  independent  replications,  (and 
equivalently  the  method  of  batch  means)  while  the  second  appendix  presents  the 
generation  of  confidence  intervals  using  the  regenerative  method  of  simulation. 

2.3  External  Control  Variates 

2.3.1  Burt,  Cover,  and  Perlas  (1970).  The  effectiveness  of  several  variance 
reduction  techniques  on  project  graph  analysis  (PERT,  GERT,  CPM,  etc.)  network 
simulations  is  examined  in  this  article.  The  techniques  considered  are  antithetic  variates, 
stratification,  and  what  we  would  now  call  external  control  variates.  They  generate  the 
external  control  variates  by  first  constructing  simplified  networks  that  have  analytical 
solutions  and  are  “similar^’  to  the  networks  under  study.  They  then  simulate  both  models 
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separately  using  the  same  random  number  stream  to  drive  both  simulations.  The 
generated  control  variates  are  then  applied  in  two  ways.  First,  no  efifort  is  made  to 
determine  an  optimal  control  variate  coefficient,  and  therefore  “-1”  is  used  to  calculate  the 
controlled  response.  They  call  this  straight  control  estimates.  Burt  et.al.  then  use  the 
regression  technique  to  estimate  the  optimal  coefficient.  Two  different  networks  are 
simulated  with  control  variates  and  results  are  reported  of  significant  variance  reduction 
with  the  regression  method  being  the  most  effective. 

2.3.2  Gaver  and  Shedler  (1971).  The  application  of  external  control  variates  to 
simulations  of  a  multiprogrammed  computer  system  is  examined  by  the  authors.  In  a 
similar  manner  to  the  study  above,  Gaver  and  Shedler  propose  a  model  of  a 
multiprogrammed  computer  system  and  a  similar  model  that  can  be  solved  analytically.  In 
the  same  marmer  as  above,  both  models  are  simulated  and  the  generated  control  variates 
are  applied  using  both  the  straight  control  variate  and  regression  to  estimate  an  optimum 
coefficient.  Results  indicate  significant  reduction  in  variance  with  the  best  results  fi’om  the 
use  of  regression. 

2.4  Internal  Control  Variates 

2.4.1.  Lavenberg,  Moeller,  and  Welch  (1982).  The  authors  propose  three  different 
internal  control  variates  (service  time  variable,  flow  variable,  and  work  variable)  that  can 
be  used  in  simulation  studies  of  a  wide  class  of  closed  queueing  networks.  The  closed 
networks  under  study  consist  of  a  finite  number  of  service  centers  with  one  or  multiple 
servers.  Additionally,  there  is  a  fixed  number  of  customers  of  different  classes.  Service 
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time  variables  are  the  sample  mean  service  time  for  a  particular  class  of  customers  at  a 
specific  service  center.  The  fi-action  of  service  completions  at  a  specific  service  center  that 
are  for  a  particular  class  of  customers  is  defined  as  a  flow  variable.  The  product  of  these 
two  variables,  the  sum  of  service  times  for  a  class  of  customers  at  a  service  center, 
constitute  the  work  variables.  Lavenberg  et.  al.  then  perform  experiments  on  a  simple 
model  of  an  interactive  multiprogrammed  computer  system.  This  same  model  is  used  in 
this  study  to  assess  the  performance  of  analytical  control  variates  on  closed  queueing 
networks.  Results  from  their  experiments  show  work  variables  produce  significant 
reduction  in  variance  over  several  different  experimental  configurations. 

2.4.2  Wilson  and  Pritsker  (1984).  Standardized  work  variables  are  developed  for 
queueing  systems  with  the  regenerative  property  in  this  study.  Given  a  service  process 
^j{ky. y  ^  ij,  at  service  center  k  with  known  expected  value  and  variance  of  cr^  , 
standardized  work  variables  are  defined  as 

a(.k,t) 

Z  K  <*)-/'.] A. 

y=l 

where  a(k,t)  is  the  number  of  service  times  started  at  center  k  during  the  time  period  [0,t]. 
These  standardized  work  variables  are  used  later  in  this  study  and  are  defined  for  use  with 
the  models  under  study  in  sections  3.3.1  and  3.4.1  below.  A  need  for  these  standardized 
variables  was  seen  since  the  work  variables  presented  in  Lavenberg  (1982)  have  been 
shown  to  have  an  asymptotic  variance  equal  to  zero  which  causes  the  variance  covariance 
matrix  for  a  set  of  these  controls  to  be  asymptotically  singular.  Since  an  inverse  of  this 
matrix  (or  usually  an  estimate  of  it)  must  be  computed  to  find  (or  estimate)  the  optimal 
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control  coefficient,  numerical  problems  can  arise  as  replication  length  increases.  The 
standardized  work  variables,  on  the  other  hand,  have  an  asymptotic  variance  of  1. 
Experiments  on  a  simple  network  were  performed  with  substantial  variance  reduction 
reported. 


2.4.3  Bauer  and  Wilson  (1993).  Another  set  of  standardized  control  variables, 
standardized  routing  variables,  are  developed  which  can  be  used  for  discrete-event 
simulation  models  that  use  a  multinomial  construct.  These  control  variates  attempt  to 
exploit  the  correlation  between  departures  from  the  mean  branching  probabilities  in  a 
network  and  the  resultant  network  response.  Standardized  routing  variables  are  defined  in 
the  following  manner.  Consider  a  multinomial  branching  process  of  g  branches  leading  to 
g  service  centers,  and  define  an  indicator  variable  as 


W)  = 


[l  if  the  z  -  th  departing  customer  goes  to  center  j, 
1 0  otherwise. 


Then  a  standardized  routing  variable  for  center  j  is  defined  as 

W)-pU) 


N{t) 


S  {Ar(/)[l-p(y)];>0)) 


1/2 


,  J=l,...,g. 


where  N(t)  is  the  total  number  of  transits  through  all  g  branches  in  the  time  interval  [0,Z] 
and  p(j)  is  the  probability  that  a  customer  is  sent  to  center  j.  Standardized  routing 
variables  are  used  in  this  study  and  are  defined  in  terms  of  the  study  model  in  section  3.4.1 
below.  Bauer  and  Wilson  performed  several  experiments  on  a  simple  network  with  results 
indicating  significant  confidence  interval  length  reduction,  particularly  when  the 
standardized  routing  variables  are  used  in  conjunction  with  standardized  work  variables. 
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2.5  Analytical  Control  Variates 

2.5.1  Tomick  (1988)  and  Tomick,  Litko,  and  Bauer  (1989).  Both  of  these  works 
report  variance  reduction  achieved  for  external  control  variates  that  are  obtained  by 
solving  approximate  queueing  models.  These  types  of  control  variates  are  called  analytical 
in  this  study  to  differentiate  them  from  external  control  variates  obtained  via  the  simulation 
of  a  separate  simulation  model.  Throughout  the  rest  of  this  discussion,  Tomick’s  external 
control  variates  will  be  called  analytical  control  variates  to  maintain  consistency  in  this 
study.  Tomick  computes  the  analytical  control  variates  two  ways.  The  first  method 
consists  of  solving  a  related  Jackson  network  using  the  observed  values  of  the  input 
random  variables  for  each  replication.  The  second  method  uses  the  Queueing  Network 
Analyzer  software  developed  by  Whitt  [22]  with  inputs  of  the  observed  values  of  the 
random  variables  for  each  replication.  Experiments  were  conducted  on  a  simple  open 
queueing  network  in  Tomick  (1988)  that  indicated  analytical  control  variates  achieved 
significant  levels  of  variance  reduction,  however  confidence  interval  coverages  were  less 
than  satisfactory  indicating  a  bias  problem  with  the  analytical  control  variates  under  study. 

2.5.2  Dietz  and  Harmonosky  (1989).  This  paper  examines  the  applicability  of  a  single 
internal  control  variate  that  is  formulated  by  aggregating  many  of  the  key  input  random 
variables  for  use  with  a  simulation  model  of  an  aircraft  sortie  generation  model.  The 
control  variate  used  is  an  analytical  model  of  the  “hands-on”  time  between  each  sortie, 
which  includes  all  ground  activities  and  flying  time  with  inputs  being  the  observed 
parameters  from  the  simulation  run.  Obviously  this  single  internal  control  variate  is  the 
same  as  the  analytical  control  variate  referred  to  in  this  study.  Several  experiments  were 
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run  for  different  values  of  aircraft  reliability,  aircraft  maintainability  and  maintenance 
resource  availability.  Results  indicate  significant  variance  and  confidence  interval  length 
reduction  can  be  obtained  using  this  approach.  Results  on  confidence  interval  coverage 
are  not  reported. 
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3.  Methodology 


3.1  Introduction 

The  performance  of  analytical  control  variates  will  be  compared  to  that  of  various 
internal  and  external  control  variates.  This  will  be  accomplished  on  several  configurations 
of  two  queuing  network  models.  All  control  variates  will  be  compared  using  the  Method 
of  Independent  Replications  [12]  as  described  in  section  3.2.  The  method  is  commonly 
referred  to  as  the  “regression  method”  since  linear  regression  formulas  can  be  used  to  find 
controlled  responses  as  well  as  estimate  controlled  variance  values  and  their  associated 
confidence  intervals. 

The  folloAving  sections  present  the  two  queueing  networks  and  the  particular 
statistical  performance  measurements  for  each  model.  The  first  section  addresses  a  simple 
open  queueing  network  and  the  second  deals  with  a  more  complex  closed  queuing 
network.  Each  section  includes  model-specific  discussion  on  the  selection  of  internal 
control  variates,  external  models  for  external  control  variates,  and  analytical  control 
variates. 

3. 2  Control  Variate  Method  of  Independent  Replications. 

The  following  methods  will  be  used  to  calculate  controlled  responses  and  estimate 
controlled  response  variances  and  confidence  intervals.  First  define  p  as  the  expected 
value  of  a  system  response  measurement  which  we  wish  to  estimate  and  Yj  as  an  unbiased 
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estimator  of  jj.  for  the  j-th  replication  of  a  simulation.  Then  define 


f,  =  r(n)=l-±Y, 


nj=, 


(3.1) 


Let  C  =  (C,,C2,...,C^)'  be  a  qxl  vector  of  control  variates  with  the  known  mean 
vector  of  ■  Let  B  =  be  a  qx  \  vector  of  constants. 

Then  a  controlled  estimate  of  an  expected  response,  /i,  is  given  by  [13] 


(3.2) 


The  variance  of  Y{B)  is  minimized  by  the  optimal  vector  of  control  coefficients  given  by 


Ecc  (3.3) 

where  is  the  covariance  vector  between  Y  and  C,  and  the  covariance  matrix 

of  C.  Since  these  values  are  normally  unknown,  p  is  estimated  by 

p^SycScc  (3-4) 

where 


Scc=^±(c,-<^c,-cy 

Then  a  controlled  estimate  of  ju  is  given  by 

Y(P)  =  nn)~P'(C-Mj 

If  we  assume  both  Y  and  C  have  the  following  joint  normal  distribution 


T 

fJ-Y 

(JjYY 

C 

i+g 

../“cl 

CY 

Lcc__ 

(3.5) 


(3.6) 


(3.7) 
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then  Y{P)  is  an  unbiased  estimator  of  jn  and  the  estimated  variance  is  given  by  [2] 


where 


and 


n  («  - 1) 


(3.8) 


(3.9) 


-SycScXc] 

n-q-l 


(3.10) 


Then  a  100(1-  a)%  confidence  interval  is  given  by 

Y{h±t,-.,t.^,-,DSy,c  (3.11) 

Since  we  have  assumed  Y  and  C  to  have  a  multivariate  normal  distribution,  the 
application  of  control  variates  using  p  can  be  considered  a  classical  regression  problem 
[9].  Computationally,  this  can  be  a  simpler  method  for  finding  Y{fi)snd.  estimating  it’s 
variance  and  associated  confidence  interval.  Since 

E[Y\C^c\  =  fx  +y3'(c-//,)  (3.12) 


the  problem  can  be  stated  as 


Y  =  Xy  +  s 

where  T^%Y„...X] 


and  Y 


P 


(3.13) 
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with 


X 


1  •  ^cx-^^c 


1  ^Xn-f^X  - 


(3.14) 


Then  to  find  the  least  squares  estimators  set 


^  ={X'XY'X’Y 


(3.15) 


Since  [11] 


Yin)  =  U  +  nC-^i^) 


(3.16) 


it  follows  that  Ji  is  equal  to  Y{P) .  Then  from  regression  theory  results  for  jii  [16],  the 


estimated  variance  of  Y{P)  is  given  by 


Var{Y{p))  =  s,,MSE 


(3.17) 


where  5n  is  the  upper  left  element  of  (XX)  and  MSE  can  be  computed  as 


MSE  = 


^  1  ^ 
\n-q-lj 


YT-rX’Y 


(3.18) 


Then  the  100(l-a)%  confidence  interval  can  be  given  by 


(3.19) 


3.3  Open  Queuing  Network. 

The  open  queuing  system  studied  consists  of  three  server  stations  connected  in 
series,  (see  figure  3.1).  Call  this  model  Mi.  Each  station,  or  node,  has  a  dedicated  single 
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server  that  can  provide  service  to  only  one  entity  at  a  time.  All  arrivals  enter  the  system  at 
node  1  with  balking  not  considered.  Upon  completion  of  service  at  node  1  all  entities 
proceed  on  to  nodes  2  and  3  consecutively.  Entities  depart  the  system  after  receiving 
service  at  node  3. 


Customer 

Arrivals  -  E^qxmential 


Service 

Center 

1 

(Infinite  queue  capacity) 


Service  time  ~  WeibuU 


Service 

Center 

2 

(Infinite  quaie  capaciTy) 

^3E0“ 

Service  time  -  WdbuU 


Service 

Center 

3 

(Queue  capacity  =2) 

Customers  leave  system 


Service  time  WeibuU 


Figure  3.1.  Open  queueing  system  model  M\ 

System  interarrival  times  are  HD  exponential  random  variables  with  a  mean 
interarrival  time  equal  to  X  (Poisson  arrival  process)  and  probability  density  function 


fix)^ 


ifx>0 


(3.20) 


0 


otherwise 


Each  service  node  has  a  different  service  time  distribution,  but  each  is  an  HD  WeibuU 
random  variable  with  probability  density  function 


/(^)  = 


^e  if  x>  0;a,P  >  0 
0  otherwise 


(3.21) 


Node  I’s  service  time  distribution  has  shape  parameter  a\  and  scale  parameter  node 
2’s  parameters  are  Oi  and  ^  ,  and  node  3’s  parameters  are  as  and  .  Then  the  mean 
service  time  for  each  node  is  equal  to 
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i=  1,2,3  (3.22) 


y  cck  J 


with  a  variance  of 


/ 

r 

1  -  r 

(a,  +0 

2\ 

V 

1  a,  J 

1 

y  a,  j 

J 

1,2,3 


(3.23) 


(Typical  service  time  and  arrival  probability  density  functions  are  shown  in  figure  3.2.) 
Nodes  1  and  2  have  an  infinite  queue  capacity  while  node  3  has  a  finite  queuing  capacity 
of  2.  If  node  3’s  queue  is  full  and  an  entity  has  finished  service  at  node  2,  that  entity  will 
wait  at  node  2  and  block  node  2’s  server  until  a  spot  becomes  available  in  node  3’s  queue. 


Figure  3.2.  Typical  service  time  distributions  (Weibull:  fx  =  10.0,  15.0,  and  18.0)  and 
typical  interarrival  distribution  (exponential:  A  =  20.0). 

This  network  is  a  simple  model  of  a  hypothetical  fast-food  restaurant  drive  through 

window  system.  The  first  node  represents  the  order  microphone,  the  second  node 
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represents  the  pay  window,  and  the  third  node  represents  the  pick-up  Avindow.  Entities 
being  served  within  the  system  are  fast-food  drive  through  customers.  The  model  is  not 
meant  to  represent  any  specific  facility,  but  has  been  created  for  illustrative  purposes  and 
the  system  assumptions  have  been  made  for  ease  of  analysis. 

The  system  performance  measurement  that  will  be  estimated  is  the  expected  time  a 
customer  spends  in  the  system  (sojourn  time)  under  steady  state  conditions.  Let  //  be  the 
true  steady  state  expected  sojourn  time  which  we  wish  to  estimate  and  let  be  the 
sojourn  time  of  the  /-th  customer  during  the y-th  replication.  Also  let  /  be  the  number  of 
customers  in  the  warm-up  period,  m  be  the  total  number  of  customers  for  each  replication, 
and  n  the  number  of  replications.  Then,  let 

1  ” 

y=l,2,,.,.n  (3,24) 

W  I  ,=;+! 

SO  that  jj.  can  be  estimated  using  the  replication/deletion  method  [13]  by 

p  =  (3.25) 

which  has  an  estimated  variance  of 


Y^(Y,-Yin)y 


n 


n 


n-  1 


and  an  approximate  100(1 -or)  percent  confidence  interval  given  by 


Y(n)±t 


n-\^l-af2 


(3.26) 


(3.27) 
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3.3.1  Internal  Control  Variates  (Open  Queueing  Network).  Internal  control  variates 
are  random  variables  generated  within  the  simulation  that  have  known  means.  For  the 
open  queueing  model,  candidate  random  variables  can  be  derived  from  the  interarrival 
distribution  and  the  three  service  time  distributions.  Standardized  work  variables,  which 
are  a  function  of  the  service  time  distribution  at  a  specific  service  center,  have  been  shown 
to  be  good  control  variate  candidates  [23].  Let  the  service  time  process  at  service  center  k 
be  the  IID  sequence  U^  (k):i  >  1,  =  1,  2,  3;  (1  ==  order  microphone,  2  =  pay  window,  3  = 
pickup  window).  Then  a  standardized  work  variable  is  defined  by 


1  "" 

=y 

n-l±, 


UAk)-ju, 


J=  1,2, ...,« 


(3.28) 


where  pk  and  ctk  are  given  by  equations  (3.22)  and  (3.23)  respectively  for  model  Mi.  Each 
Wkj  has  a  mean  of  zero  as  (/w  -  /)  — >  oo  [23]. 

Another  candidate  internal  control  variate  is  a  standardized  interarrivaJ  variable. 
Similarly  let  ;  /  >  1,  be  the  IID  sequence  of  system  interarrival  times.  Since  X  is  the 
mean  of  the  interarrival  time  distribution  for  model  Mi,  define  a  standardized  interarrival 
variable  as 


^  ■yim-l  i. 


V  -X 


-i+i 


cr„ 


1,2,  ...,n 


(3.29) 


which  also  has  a  mean  of  zero  as 

Given  the  four  candidate  internal  control  variates  defined  above,  fifteen  (2^  -  1) 
different  combinations  of  control  variates  could  be  used  to  reduce  the  variance  of  the 
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estimated  average  sojourn  time.  Let  be  the  matrix  of  all  candidate  internal  control 
variates  such  that 


Wn  .. 

^2.  - 

fV 

W,,  .. 

■■ 

A,  .. 

..  A„ 

Also  let  and  A  =  .  Then  let 


j=t 


(3.30) 


fFi 

W2 

Ws 

A 


(3.31) 


be  the  vector  of  internal  control  variate  means.  Then  for/?  =  1,2,  ...,  15,  is  one  of 

— JNT 

the  possible  matrices  of  selected  control  variates  and  Cp  is  the  associated  vector  of 
sample  means.  Table  1.1  enumerates  the  different  internal  control  variate  combinations. 

A 

Let  Y(P)^  be  the  internally  controlled,  estimated  expected  sojourn  time  using  ,p  = 
1, 2, ...,  15,  so  that 

^  ^  ^^INT 

Y(py^  =Y(n)-P'(Cp  -0)  (3.32) 

Then  each  Y(P)y^  can  be  calculated  using  equation  (3.15)  and  the  estimated  variance 

and  confidence  intervals  can  be  found  using  equations  (3.17)  and  (3.19)  respectively. 
Performance  of  each  controlled  estimator  can  be  compared  based  on  variance  reduction 
and  confidence  interval  reduction  for  each  estimator. 
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p 

w.vr 

L,  p 

1 

\WA 

2 

[W2\ 

0 

{WA 

4 

\A\ 

5 

[Wx  W2V 

6 

\WxW,V 

7 

\WxAV 

8 

[W2WA' 

9 

{W2AV 

10 

{W^AV 

11 

[Wi  W2W3V 

12 

[  Wi  fV2Ay 

13 

1  Wx  WsAV 

14 

1  W2W3Ay 

15 

[  Wx  W2  WiAy 

Table  3.1.  All  possible  combinations  of  internal  control  variates 
3.3.2  Analytical  Control  Variate  (Open  Queueing  Network).  Application  of 
anal5^ical  control  variates  requires  the  development  of  a  separate  model,  say  M2,  that  can 
be  solved  analytically.  Additionally,  the  solution  to  the  analytical  model  M2  must  be 
correlated  to  the  response  of  each  replication  of  the  simulation  model  Mi  in  order  to 
reduce  the  variance  of  the  response  using  control  variate  techniques.  In  this  specific  case, 
the  analytical  model  is  a  series  M/M/l  open  queuing  network  with  three  nodes  as  depicted 
in  figure  3.3.  The  interarrival  rate  is  MX,  same  as  Mi,  and  the  three  service  rates  are  given 
by  I///1 ,  \lp2 ,  and  I///3,  with  pk  given  by  equation  (3.21).  Based  on  parameters  a*  and pk 
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from  model  M\.  Then,  from  queuing  theory,  the  expected  steady  state  sojourn  time  is 
known  to  be  [18] 


.ANALTT 


1 


J_ 


+  - 


1 


_L_1 

//2  X 


+ 


1 


//3 


J 


(3.33) 
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Figure  3.3.  Open  M/M/ 1  queueing  network  model  M2  and  M3 
Define  as  the  1 X  M  vector  of  analytic  control  variates  for  the  n  replications. 

To  find  ,  first  define  Tij  as  the  /-th  customer’s  interarrival  time  for  the  y-th 

replication  of  Mi.  Using  the  replication/deletion  method  with  a  warm-up  of  /  customers 
and  a  total  of  m  customers  simulated,  estimate  the  expected  interarrival  time  during  the  j- 
th  replication  by 

(3.34) 

^  / ;=/+] 

Similarly  define  Vij{k)  as  the  service  time  for  the  /-th  customer  at  the  k-t\\  service  center 
during  the  y-th  replication  of  Mi.  Then  to  estimate  the  expected  service  time  at  the  ^-th 
service  center  during  the  y-th  replication  let 

V j{k)^-l-  ^v  {k)  1,2,  3;  y=  1,2,  ...,//  (3.35) 


3-11 


Then  ,  the  analytic  control  variate  for  the y-th  replication  is  calculated  as 


with 


C 


ANALTT 


1  1 


-  +  - 


1 


1 


1 


Vji\)  Rj  Vj(2)  Rj  Fj(3)  Rj 


(3.36) 


-ANALYT 


ANALl'T 


(3.37) 


Define  Y{P)  as  the  analytically  controlled  estimate  of  expected  sojourn  time. 


where 


Yf^p-^ANAUT  ^  Y{n)-p\C  -^“^) 


(3.38) 


which  can  be  computed  using  equation  (3.15).  The  estimated  variance  and  confidence 


interval  of  Y{P)  can  be  calculated  using  equations  (3.17)  and  (3.19). 


3.3.3  External  Control  Variate  (Open  Queueing  Network).  Application  of  external 
control  variates  requires  a  separate,  similar  simulation  model  with  a  known  expected 
sojourn  time,  Call  the  external  simulation  model  M3.  Using  common  random 

numbers  between  models  M\  and  M3  the  estimated  expected  sojourn  time  between  the  two 
models  should  be  highly  correlated  so  that  variance  reduction  can  be  obtained  using 
control  variate  techniques. 

For  our  example,  simulation  model  M3  has  the  same  characteristics  and  parameters  as 
model  M2  —  a  series  of  three  M/M/1  infinite  capacity  queues  with  the  same  interarrival  and 
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service  time  distributions  as  M2  (see  figure  1 .2).  Therefore 


.EXJ 


.ANALYT 


1 


Ml  ^ 


■  +  ■ 


1 


J_ 

Ml 


+ 


1 


Ml  ^ 


(3.39) 


Random  numbers  are  synchronized  between  Mi  and  Ms  for  all  arrival  times  and 
service  times  at  each  service  center.  Since  both  Mi  and  M2  have  exponential  interarrival 
times,  customer  arrival  times  will  be  the  same  for  each  simulation  model.  At  the  service 
centers,  the  same  random  number  streams  of  uniformly  distributed  random  variables 
between  zero  and  one  are  used  to  generate  the  random  service  time.  The  streams  will 
produce  different  service  times  since  the  distributions  are  different  for  each  model 
(Weibull  and  exponential  respectively),  but  will  be  highly  correlated  since  both  Weibull 
and  exponential  random  variables  are  generated  using  the  inverse  transform  method  [17] 
and  therefore  have  a  one  to  one  relationship  with  each  uniform  [0,  1]  random  variable 
generated  by  each  stream. 

To  generate  the  control  variates,  let  yf^  be  the  system  sojourn  time  for  the  /-th 

customer  during  the y-th  replication  of  external  model  Ms.  As  with  model  Mi,  m  customer 
arrivals  are  simulated  (with  a  warm-up  period  of  I  customers),  in  order  to  find  the  mean 
sojourn  time  for  replication  j  as 

yT=~±yr  (340) 

m-ltt 


Then  define  Cf^,  the  external  control  variate  for  the  y-th  replication,  as 


=  Yj 


—ext 


j  =  1,  ...,n. 


(3.41) 
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and 


f^EXT 

Then,  the  externally  controlled  estimated  sojourn  time  is  given  by 

yCP)^  =  Y{n)  -  ) 


c""=-T 

■MXt 


n 


/=! 


(3.42) 


(3.43) 


which  can  be  computed  using  equation  (3.15).  As  with  the  other  controlled  responses,  the 
estimated  variance  and  confidence  interval  can  be  found  using  equations  (3.17)  and  (3.19). 


3.4  Closed  Queueing  Network 

The  closed  network  studied  has  the  basic  form  shown  in  figure  3.4.  The  network 
structure  is  the  same  as  the  network  examined  by  Lavenberg,  Welch,  and  Moeller  (1982). 
An  additional  variation  of  this  network,  (shown  in  figure  3.5)  will  be  examined. 


Figure  3.4.  Closed  queueing  network  model  Qi 
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As  depicted  in  figure  3.4,  the  network  has  a  total  of  S  service  centers  Avith  a  total  of 
TV  customers  within  the  closed  network.  Call  this  model  Qi.  Service  center  1  has  exactly 
N  servers  resulting  in  no  queues  at  center  1 .  The  remainder  of  the  centers  are  single  server 
queues  with  customers  being  served  with  a  ‘first  in  -  first  out’  queue  priority.  The 
underlying  transition  probability  matrix  for  customer  movement  throughout  Qi  is  given  by 

'O  1  0  ..  o' 

Pi  0  ••  Ps 

0  1  0  ..  0  (3.44) 

0  1  0  ..  0 

with  >0,  k  ^  1,  3,  ...,  S.  As  indicated  above,  all  customers  completing  service  at 
center  1  are  routed  to  center  2.  Upon  completion  of  service  at  center  2,  the  customer  is 
routed  to  centers  1,  3, ...,  or  S  with  probability  pk,k=\,3, 

The  service  time  distributions  at  centers  1,  3,  ...,  S  are  characterized  as  different  IID 
Weibull  random  variables  with  probability  density  function 


fix)  =  P 

0 

V 


if  x>  0-,a,P  >  0 
otherwise 


(3.45) 


Each  of  these  centers  has  a  service  time  distribution  with  shape  parameter  a*  and  scale 
parameter  Pk,k=  \,3,  ..,S.  Then  the  mean  service  time  for  each  service  center  is  equal  to 

(3.46) 

\  a,  J 
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and  with  variance  given  by 


=  Pt 


V  a.  J 


-r 


y  CC,  J 


2^ 


k=l,3,...,S  (3.47) 


Service  center  2  has  an  exponential  service  time  distribution  with  probability  density 
function 


/(^)  = 


—e 

A 

0 


ifx>0 

otherwise 


(3.48) 


with  mean  service  time  X. 

Q\  can  be  considered  a  simple  model  of  an  interactive  multiprogrammed  computer 
system  [10].  The  customers  in  this  network  are  users  of  the  system  with  service  center  1 
representing  the  user  terminals  with  the  a  service  time  associated  with  the  user’s 
“thinking”  time  between  system  requests.  Service  center  2  represents  the  system’s  central 
processing  unit  (CPU)  and  service  centers  3,  ...,S  model  the  system’s  mass  storage  units 
(disk,  drum,  tape,  etc  ).  The  service  time  at  center  2  represents  processing  time  for  the 
user’s  task  request  until  either  the  task  is  completed,  which  will  cause  the  customer  to 
return  to  center  1  and  start  another  “think”  time,  or  until  data  from  a  mass  storage  unit  is 
required.  Service  times  at  centers  3,  ...,S  model  the  time  required  to  transfer  data  from 
the  device  it  represents  to  the  main  memory  where  it  can  be  acted  upon  by  the  CPU 
(center  2). 


3-16 


Since  users  must  be  allocated  a  portion  of  main  memory  in  order  to  access  the  CPU 
and  storage  devices  in  these  types  of  systems,  main  memory  limitations  may  not  make  it 
possible  for  all  users  to  have  access  at  the  same  time.  This  leads  to  the  variation  of  Q\ 
discussed  previously,  called  Qc,  as  shown  in  figure  3.5.  In  Oc  at  most  N'  <  A'" customers 
can  enter  the  subnetwork  of  centers  3,  ...,  S  with  service  center  3  representing  the  CPU 
and  centers  4, ...,  5”  being  the  mass  storage  devices.  Service  center  2  holds  customers  until 
the  number  of  customers  in  the  subnetwork  is  less  than  N' .  The  service  time  random 
variables  for  center  2  are  identically  zero  and  there  is  no  queueing  at  center  2  if  the 
number  of  customers  in  the  subnetwork  is  less  than  N’ . 


r 


Figure  3.5.  Closed  queuing  network  model  Qc 
Several  measurements  of  system  performance  may  be  of  interest  for  such  a  system, 
but  we  will  focus  on  two  of  them.  The  first  being  the  steady  state  expected  system 
response  time,  defined  as  the  long  run  average  time  it  takes  for  customers  departing 
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service  center  1  to  return  to  center  1 .  The  other  performance  measurement  of  interest  is 
the  steady  state  expected  CPU  utilization  —  the  fraction  of  time  service  center  t\vo  is 
serving  customers.  It  should  be  noted  that  these  performance  measurements  are  of 
interest  in  both  Oi  and  Oc  and  the  following  discussion  on  estimating  these  values  applies 
to  both  models. 

To  estimate  the  expected  system  response  time  for  Oi,  using  computer  simulations, 
let  r  be  the  true  steady  state  expected  system  response  time  which  we  wish  to  estimate, 
and  let  ty  be  the  response  time  of  the  /-th  return  to  center  1  during  the  y-th  simulation 
replication.  If  we  define  an  event  as  the  completion  of  service  at  any  of  the  service 
centers,  terminate  each  replication  upon  completion  of  M  events.  Let  Z  be  the  number  of 
events  in  a  warm-up  period  in  order  to  eliminate  bias  from  the  initial  transient  behavior. 
Define  nij  as  the  number  of  times  customers  return  to  center  1  before  event  M  and  Ij  as  the 
first  return  to  center  1  following  event  Z  for  replication  j.  Then  let 


_J__V 


7  =  1,2, . ,n 


(3.49) 


so  that  rcan  be  estimated  using  the  replication/deletion  method  [13]  by 


?  =  r(«)  =  -2^ 


(3.50) 


which  has  an  estimated  variance  of 


TiT^-T{n)f 

1 


n  n  n  —  \ 


(3.51) 
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and  an  approximate  100(1 -or)  percent  confidence  interval  given  by 


7’(w)  ±  i_„/2  ■ 


(3.52) 


CPU  utilization  can  be  estimated  in  the  following  manner.  Let  v  be  the  true  steady 
state  expected  CPU  utilization.  Now  let  by  be  the  amount  of  simulated  time  service  center 
2  is  serving  the  /-th  customer  to  visit  center  2  during  replication  j.  Define  ej  as  the 
simulated  time  of  replication  j  when  event  M  occurs  and  gj  as  the  simulated  time  of 
replication  j  when  event  Z  occurs.  Further  let  qj  equal  the  number  of  customers  that 
complete  service  at  center  2  before  time  Cj  and  let  rj  be  the  first  customer  to  complete 
service  at  center  2  following  time  gj  during  replication  j.  Then  let 


1 


j  =  l,2,...,n 


(3.53) 


so  that  V  can  be  estimated  using  the  replication/deletion  method  [13]  by 


(3.54) 


which  has  an  estimated  variance  of 


n  n  n-\ 


(3.55) 


and  an  approximate  100(1 -a)  percent  confidence  interval  given  by 


U{ri)  ±  ■ 


(3.56) 
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3.4.1  Internal  Control  Variates  (Closed  Queueing  Network).  Possible  random 
variables  with  known  means  for  model  0\  include  the  S  service  time  distributions  and  5-1 
values  of  the  routing  variables ....,ps.  Again,  the  following  discussion  applies  to  Qc 
as  well  as  0\.  As  in  the  open  queueing  network,  standardized  work  variables  should  be 
good  control  variate  candidates.  Another  set  of  candidate  internal  control  variates  with 
demonstrated  utility  that  are  appropriate  to  0\  are  standardized  routing  variables  [2], 
Each  of  these  internal  control  variates  are  applied  to  estimate  both  expected  response  time 
and  expected  CPU  utilization. 

The  standardized  work  variables  are  derived  in  the  same  manner  as  they  were  for 
model  Ml.  Let  s^ik)'.  /  >  1,  k=  1,2,  ....,  S  represent  the  service  time  process  at  the  A:-th 
service  center.  Further  let  aj{k)  be  the  number  of  service  times  completed  at  center  k,  and 
let  h/k)  be  the  first  service  time  completed  at  center  k  after  time  gj  for  replication  j.  Then 
the  standardized  work  variables  are 


.Jaj  (k)  -bj  (k)  +  l  i=bji.k)  ^ 


(3.57) 


for  j  =  1,  2, ...,  n  and  k=\,2, ...,  S 

where  /4  and  Sk  are  given  in  equations  (3.46)  and  (3.47)  respectively.  Each  Wjg  has  a 
mean  of  zero  as  aj(k)  -  bj{k)  goes  to  infinity  [23]. 
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Standardized  routing  variables  are  developed  in  the  following  manner.  Given  the  five 
branching  processes  at  service  center  2,  (center  3  for  Oc),  define  an  indicator  variable  as 


4W  = 


[  1  if  the  i  -  th  departing  customer  goes  to  center  k  for  replication  j, 
0  otherwise. 


(3.58) 


Then  a  standardized  routing  variable  for  activity  k  for  replication  j  is  defined  as 

IAk)-p, 


*J<»  I"/  (2)  -  (2)  + 1  Jl  -  p,  ]pj  I  ^ 


(3-59) 


for  1,  3,  4, ...,  .S'  and  j  =  1,  2, ...,  w 

where  [a,(2)  -  bj(2)  +1]  is  the  total  number  of  service  completions  at  service  center  2  for 
replication  j  during  the  time  interval  [gj,  Cy]  and  the  pk  are  the  transition  probabilities  from 
service  center  2  to  center  k  as  given  by  Qis  transition  probability  matrix  P.  It  has  been 
shown  that  the  vector  of  routing  variables,  R,  converges  in  distribution  to  a  multivariate 
normal  distribution  vrith  a  mean  of  0  as  the  simulation  run  length  increases  [2]. 

As  with  model  Mi,  the  25-1  internal  control  variates  make  2^’^  -  1  different 
combinations  of  control  variates,  so  define  as  the  (25  -  l)xw  matrix  of  all  possible 
internal  control  variates  for  n  replications.  Then  for;?  =  1,  2, ...,  25  -  1;  let  be  one  of 

— INT 

the  possible  matrices  of  selected  control  variates  and  let  Cp  be  defined  as  the  mean 


vector  of  the  associated  matrix  (same  as  equation  (3.30)). 

The  controlled  system  response  time,  T{P)^ ,  is  defined  by 


which  can  be  calculated  using  equation  (3.15)  and  the  estimated  variance  and  confidence 
intervals  can  be  found  using  equations  (3.17)  and  3(.  1 9). 

Similarly,  the  controlled  CPU  utilization  fi-action  U (yff)*^  is 

Uihf  -  U{n)  -  P’(C7  -  0)  (3.61) 

and  the  associated  variance  and  confidence  intervals  can  be  estimated  using  the  regression 
equations  of  section  5.2. 

3.4.2  Analytical  Control  Variate  (Closed  Queueing  Network).  To  find  anal5d:ical 
control  variates  for  both  Qi  and  Qc  consider  a  new  model,  O2,  that  can  be  solved 
analytically  and  is  similar  in  structure  to  models  Qi  and  Qc-  Let  O2  have  the  same 
structure  and  transition  probability  matrix  as  model  Qi,  except  that  all  service  time 
distributions  for  O2  are  HD  exponential  random  variables  with  mean  service  times  equal  to 
those  of  Qi.  Using  the  Mean  Value  Analysis  (MV A)  algorithm  [11]  the  steady  state 
expected  system  response  time  and  CPU  utilization  fraction  can  be  determined  exactly. 
Let  and  be  the  expected  response  time  and  expected  CPU  utilization  for  the 

analytical  model  02 

The  MVA  algorithm  is  for  closed  networks  that  have  closed  form  solutions.  The 
algorithm  yields  the  mean  values  of  the  network  performance  characteristics  response 
time,  queue  length,  throughput,  and  utilization  for  each  service  center  in  the  network. 
Basically  the  MVA  solves  the  network  by  first  determining  system  characteristics  with 
only  one  customer  in  the  system  and  then  solves  it  for  two  customers  based  on  the 
information  obtained  for  one  customer  using  the  mean  value  theorem.  The  mean  value 
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theorem  relates  the  response  time  of  a  service  center  when  n  customers  are  present  to  the 
length  of  the  queue  at  the  service  center  when  « - 1  customers  are  present.  The  algorithm 
is  reapplied  until  the  system  is  solved  for  customers. 

To  calculate  the  value  of  the  analytical  control  variate  for  each  replication,  the 
ForQue  computer  program  [7],  which  implements  the  MVA  algorithm,  is  used.  ForQue 
requires  the  following  information  about  a  network  in  order  to  arrive  at  a  solution; 
number  of  customers  {N),  number  of  service  centers  (5),  number  of  servers  at  each  center 
{N  at  center  1,  one  at  all  the  others),  mean  service  times  at  each  center  (/4),  and  the 
transition  probability  matrix  (P).  Obviously,  N,  S,  and  the  number  of  servers  remain 
constant  for  each  replication.  Additionally,  rows  1,  3, ...,  5  of  P  remain  the  same  for  each 

replication  and  are  all  identically  equal  to  the  1x5"  vector  [0  10 . 0].  To  find  the 

analytical  control  variates  the  mean  service  times  and  actual  branching  proportions  from 
center  2  must  be  calculated  and  input  to  ForQue  for  each  replication. 

Define  as  the  analytical  control  variate  for  expected  response  time  for 


the  j-th  replication.  Also  let  be  the  analytical  control  variate  for  expected 

CPU  utilization.  To  find  each  of  these  control  variates  using  ForQue,  define  Vij{k)  as  the 
service  time  of  the  /-th  customer  to  visit  center  k  during  replication  j.  Then  to  find  the 
average  service  time  at  each  center  let 


mj(k) 

7 


{k)  A:=l,2,...,5';y=l,2,....,«  (3.62) 
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where  mj{k)  is  the  number  of  customers  completing  service  at  center  k  and  //A:)  is  the  first 
customer  to  complete  service  at  center  k  following  time  gj  for  replication  j.  To  find  the 
realized  routing  proportions  fi-om  center  2,  define  the  indicator  variable  h/k)  that  equals 
one  when  the  /-th  customer  leaving  center  2  enters  center  k,  and  zero  otherwise  for 
replication  Then  let 

mj{k) 

Z4(*) 

Pf  =  a  -  *=1,3,  ..,5  andi-1,2,  ...K  (3.63) 

S  Z4W 

t=l  i=Ij(2) 

Given  these  values,  following  each  replication  J,  input  Vjk  for  k=  \,2, S  and  pj^  for  k 
=  1,  3,  S  into  ForQue  to  obtain  and  .  Then  the  controlled 


responses  are 


and 


jf^P^ANALYT  ^ 

JJ^PyANALYT  ^  _  p,^^jj^ANALYT  _  ^ANALYT  ^ 


(3.64) 

(3.65) 


where  and  are  the  means  of  the  n  values  of  each  analytical 


control  variate.  These  controlled  responses  and  their  associated  estimated  variances  and 
confidence  intervals  can  be  calculated  using  the  equations  of  section  3.2. 

3.4.3  External  Control  Variate  (Closed  Queueing  Network).  External  control 
variates  can  be  obtained  by  creating  a  simulation  model,  O3,  of  the  same  form  as  model 
O2.  The  true  steady  state  expected  response  time,  and  expected  CPU  utilization, 
,  will  be  the  same  as  those  for  model  O2.  Using  common  random  numbers,  this 
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simulation  model  can  then  be  used  to  produce  external  control  variates  for  both  models  Qi 
and  Qc.  However,  due  to  the  nature  of  a  closed  network,  exact  synchronization  of 
common  random  numbers  becomes  very  difficult. 

Since  the  service  time  distributions  of  centers  1,  3, ...,  5'  in  Oi  are  IH)  Weibull  random 
variables,  it  is  possible  to  use  the  same  underlying  uniform  random  number  to  produce  the 
/■-th  service  time  at  service  center  k  in  O3,  producing  highly  correlated  service  times. 
Additionally,  the  routing  random  variables  for  the  /-th  customer  to  complete  service  at 
service  center  2  can  be  the  same  random  number  for  both  models  0\  and  Qi-  The  inability 
to  produce  exact  synchronization  with  this  scheme  comes  from  the  fact  that  the  /-th 
service  time  at  service  center  k,  although  highly  correlated,  will  not  be  exactly  the  same  for 
both  models  which  results  in  a  different  sequence  of  specific  customers  arriving  at  service 
center  2  for  each  model  so  that  the  each  customer  will  likely  have  a  different  sequence  of 
service  centers  visited  for  each  model.  It  seems  the  only  way  around  this  is  to  assign  all 
routing  and  service  time  random  variables  to  each  specific  customer  at  the  beginning  of 
each  replication  for  both  Qi  and  O5.  Although  this  is  possible,  it  is  an  extremely 
cumbersome  way  to  program  both  simulation  models. 

Although  exact  synchronization  of  common  random  numbers  is  not  achieved,  the 
following  scheme  for  random  number  generation  is  used  for  O3.  The  same  random 
number  streams  of  uniformly  distributed  random  numbers  between  zero  and  one  will  be 
used  to  generate  the  /-th  service  time  exponential  random  variable  for  O3  and  Weibull 
random  variable  for  Qi,  at  service  centers  1,  3,  ...,  S.  The  stream  of  routing  random 
variables  will  be  identical  for  both  models  as  well  as  the  stream  of  random  service  times  at 
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service  center  2.  Based  on  this  scheme,  the  output  of  O3,  should  still  have  some 
correlation  with  Oi  and  produce  variance  reduction. 

To  generate  the  external  control  variate  for  expected  response  time,  let  be  the 

response  time  for  the  i-th  customer  for  replication  J  of  model  O3.  Using  the  same  number 
of  events  (M)  to  terminate  O3  and  Oi,  and  the  same  number  of  events  (Z)  to  define  the 
warm-up  period  for  both  models,  define  ef^  as  the  simulated  time  in  which  M  events 

occur  for  replication  j  of  O3.  Further  define  as  the  simulated  time  Z  events  occur  for 
replication  j  of  Q3.  Then  define  as  the  number  of  customers  that  return  to  service 
center  1  before  time  and  let  be  the  first  customer  to  return  to  center  1  following 
time  gf^  for  replication  j.  Then  define  the  external  control  variate  for  response  time  as 


. 1.2, 

/W  —  I J  +1  j-fEXT 


(3.66) 


For  the  external  control  variate  for  expected  CPU  utilization,  let  be  the  amount 
of  time  it  takes  to  service  customer  i  at  service  center  2  for  replication  j  of  model  Q3. 
Also  let  qf^  be  the  number  of  customers  that  complete  service  at  service  center  2  before 

ef^  and  let  be  the  first  customer  to  complete  service  at  center  2  following  time 
g^  for  replication  j.  Then  the  external  control  variate  for  CPU  utilization  for  replication 
j  is  defined  as 


„Exr 


^EXT 


-Sj 


EXT 


z*r  y=>.2. 


n 


(3.67) 
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The  externally  controlled  responses  are  given  by 

tCP)^  =  T{n)  -  p\C(T)^  -  T^)  (3.68) 

and 

U(P)^  =  U{n)  -  P’(C{U)^  -  v^)  (3.69) 

where  C(T)^  and  C(U)^  are  the  means  of  the  n  values  of  the  respective  external 

control  variates.  These  controlled  responses  and  associated  estimated  variance  and 
confidence  intervals  can  be  found  using  the  regression  formulas. 
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4.  Results 


4.1  Introduction 

This  chapter  details  the  results  of  several  experiments  on  both  the  open  and  closed 
queueing  networks  described  in  chapter  3.  First,  experimental  procedures  used  for  both 
networks  are  discussed  below,  followed  by  a  discussion  on  computer  implementation  of 
the  simulation  models  and  the  controlled  response  formulas.  Finally,  experimental  results 
from  the  from  the  two  queueing  models  are  presented. 

4.2  Experimental  Procedures 

Several  experiments  are  conducted  on  each  network.  Number  of  replications  and 
network  parameters  are  changed  for  each  design  point.  Three  different  network  settings 
have  been  chosen  for  the  open  network  and  twelve  different  network  settings  have  been 
selected  for  the  closed  network  (six  for  Qi  and  six  for  Qc)-  Selection  of  particular 
network  settings  and  their  values  are  discussed  below  in  section  4.2.2.  At  each  network 
setting,  100  experiments  are  conducted  with  replication  size  of  10,  and  50  experiments  are 
conducted  with  replication  size  set  at  20. 

For  eveiy  experimental  setting,  performance  measures  for  internal,  analytical,  and 
external  control  variates  are  compared.  Estimated  variance,  confidence  interval  size, 
coverage,  and  estimated  Mean  Square  Error  values  for  the  controlled  responses  are 
compared  to  the  same  values  for  the  uncontrolled  responses.  For  the  internal  control 
variates,  comparisons  are  made  for  all  possible  combinations  of  internal  control  variates 
for  the  closed  network.  All  possible  combinations  of  standardized  work  variables  and 
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standardized  routing  variables  are  considered  separately  for  the  open  network.  The 
internal  control  variates  are  then  screened  and  all  combinations  of  the  best  performing 
internal  control  variates  are  considered  for  the  closed  network.  Performance 
considerations  are  described  below. 

4.2.1  Comparison  Techniques.  Comparisons  were  made  using  the  generalized 
method  presented  by  Bauer  and  Wilson  [2].  Let  p  be  the  expected  value  of  the 
performance  measurement  of  concern.  For  the  n  =  10(20)  rephcations  for  the  A^-th 
experiment,  A  =  1, 2, ...,  m{m=  100(50)),  an  estimate  of  p  is  computed.  Call  the  estimate 
Pkif),  where  /  =  1  or  2.  When  /  =  1,  no  controls  are  used  and  1-2  when  some  set  of 
control  variates  are  used.  In  a  similar  manner,  let  <7^(/)be  the  uncontrolled  (/  =  1)  and 
controlled  (/  =  2)  estimates,  for  the  A-th  experiment,  of  the  variance  of  //*(/)•  Then  the 
average  variance  estimator  over  all  m  experiments  for  a  given  setting  is 


=  '=1,2 


m 


(4.1) 
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The  percentage  change  in  estimated  variance  due  to  the  use  of  control  variates  is  then 
estimated  by  10o|a^ (1)  -  (^)]/ ^^(0  • 

For  the  k-th  experiment  the  confidence  interval  estimate  is  given  by 


K,{I)  =  p,{I)±H,{r) 


(4.2) 


where  is  the  estimated  half-width  as  given  in  equation  (3.19)  vrith  a  =  0.10.  By 


letting  L).  (/)  =  2H^  (/)  be  the  estimated  width  of  confidence  interval,  we  can  find  the 
average  width  of  the  confidence  interval  estimator  over  all  m  experiments  for  a  given 
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setting  as 


^  1  m 

L{l)  =  -Y,k{I)  1=1,2.  (4.4) 

Then,  as  with  the  variance  estimates,  the  percentage  change  in  estimated  confidence 
interval  width  due  to  the  use  of  a  particular  set  of  control  variates  is  estimated  by 

ioo[i(i)-i(2)]/Z(i). 

An  area  of  concern  for  control  variate  performance  is  the  amount  of  bias  in  the 
controlled  estimates  of  Any  bias  induced  is  due  to  the  fact  that  P  must  be  estimated  and 
is  generally  not  independent  of  Yim)  [13],  One  measurement  that  is  effected  by  bias  is  an 
estimated  confidence  interval  coverage  probability.  To  find  one,  let 

4(/)  =  P  if/^eA,(/), 

[o  otherwise 

for  /  =  1, 2  and  A:  =  1,  2, ...,  m.  As  before  let  /  =  1  for  no  controls  and  1  =  2  when  some  set 
of  controls  is  used.  Then  an  estimate  of  the  confidence  interval  coverage  probability  is 
given  by  the  actual  coverage  fi'action  for  A^  (/)  given  by 


1 

/(o=-2;4<o  '-1,2 

m 


(4.4) 


Realized  coverage  may  not  always  be  the  best  indicator  of  bias.  For  example,  a  point 
estimate  may  be  “close”  to  //,  but  if  the  associated  confidence  interval  is  small  enough, 
coverage  may  not  be  realized.  In  order  to  measure  this  “closeness”,  another  measure  of 
bias,  the  estimated  value  of  the  Mean  Square  Error  of  a  point  estimator,  will  be  computed. 
To  estimate  the  average  value  of  the  Mean  Square  Error  for  all  m  experiments  of  a  given 
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setting  let 


1  ^ 

=  /=1.2  (4.5) 

4.2.2  Open  Queueing  Network  Settings.  The  open  queueing  network  presented  in 
section  3.3  is  examined  at  three  different  settings.  For  a  given  arrival  rate  of  20  customers 
per  time  period,  the  traffic  intensity  at  every  server  station  is  set  to  approximately  0.50, 
0.75,  and  0.90  for  each  experimental  setting.  Specific  system  settings  for  model  M\  are 
given  in  table  4.1.  At  each  of  these  settings,  results  are  obtained  for  replication  sizes  of  10 
and  20.  For  the  first  two  design  points  1,000  customers  are  simulated  for  every 
replication  with  data  from  the  first  300  customers  ignored  in  order  to  minimize  the  effect 
of  the  initial  transient.  Statistics  are  gathered  on  the  final  700  customers.  More  customers 
were  simulated  for  the  third  design  point  (0.90  traffic  intensity)  since  it  takes  longer  to 
exhibit  steady-state  behavior  due  to  the  saturated  state  of  the  system.  In  this  case  20,000 
customers  transit  through  the  system  with  data  ignored  on  the  first  10,000  customers. 


Design  Point 

Activity 

Distribution 

a 

P 

Mean 

Variance 

1. 

Arrivals 

~ 

20.0 

400.000 

Center  1  service 

Weibull 

3.6150 

6000.0 

10.0 

9.448 

Center  2  service 

Weibull 

3.5416 

6000.0 

10.5 

10.808 

Center  3  service 

E^iSOlHlllii 

3.2080 

3000.0 

10.0 

13.829 

2. 

Arrivals 

~ 

~ 

20.0 

400.000 

Center  1  service 

WeibuU 

3.0850 

6000.0 

15.0 

28.266 

Center  2  service 

Weibull 

3.0489 

6000.0 

15.5 

30.830 

Center  3  service 

Weibull 

2.8355 

3000.0 

15.0 

32.883 

3. 

Arrivals 

~ 

20.0 

400.000 

Center  1  service 

Weibull 

3.4398 

30000.0 

18.0 

33.479 

Center  2  service 

Weibull 

3.4081 

30000.0 

18.5 

35.958 

Center  3  service 

Weibull 

3.7790 

3000.0 

18.0 

34.591 

Table  4.1.  Experimental  Settings  for  open  queueing  network 
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4.2.3  Closed  Queueing  Network  Settings.  For  the  closed  queueing  network  presented 
in  section  3.4.,  six  different  experimental  design  points  are  selected  for  both  models  0\ 
and  Qc.  For  all  design  points,  the  number  of  service  stations,  5,  is  6  (5"  ==  7  for  model  Qc), 
and  the  number  of  customers,  N,  is  equal  to  25.  For  model  Qc,  the  number  of  customers 
allowed  into  the  subnetwork,  N’ ,  is  equal  to  5  for  all  settings.  To  create  the  6  different 
settings  for  each  model,  2  different  transition  probability  matrices  are  applied  to  3  sets  of 
service  time  distributions.  By  applying  equation  (3.43),  the  2  transition  probability 
distributions  are  provided  in  table  4.2  and  the  3  service  time  settings  are  listed  in  table  4.3. 
Then  the  12  different  design  points  for  the  closed  queueing  network  are  outlined  in  table 
4.4.  Results  are  obtained  for  each  design  point  for  replication  size  of  both  10  and  20. 

In  all  12  cases,  replications  of  the  closed  queueing  network  are  terminated  following 
the  completion  of  2,000  events.  To  minimize  the  effect  of  the  initial  transient  behavior, 
data  from  the  first  500  events  is  ignored.  In  addition,  the  initial  state  of  the  network 
(number  of  customers  in  each  service  center)  is  based  on  the  expected  number  of 
customers  at  each  center  for  analytical  model  Qj.  These  expected  values  are  determined 
by  solving  the  system  with  the  ForQue  program.  The  probability  that  a  customer  is  at  a 
given  service  center  is  determined  by  dividing  the  expected  number  of  customers  at  a 
service  center  by  the  number  of  total  customers  (25).  Then  at  the  start  of  each  replication 
(for  Q\,  Qc,  and  O3)  each  customer  is  assigned  a  uniform  random  number  in  the  interval 
[0, 1]  and  is  routed  to  a  particular  service  center  by  the  above  probabilities. 
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Table  4.2.  Transition  probability  matrix  values  for  closed  queueing  network 


Model  Qy  (Qc) 


Service  Center 

Distribution 

a 

_ P _ 

Mean 

Variance 

Setting  A. 

1 

WeibuU 

1.46824 

1000.0 

100.00 

4795.784 

~ 

~ 

1.00 

1.000 

WeibuU 

5.6476 

10.0 

1.39 

0.081 

4(5) 

WeibuU 

5.64760 

10.0 

1.39 

0.081 

waBmm 

WeibuU 

2.61249 

1000.0 

12.50 

26.442 

WeibuU 

2.61249 

1000.0 

12.5 

26.442 

Setting  B. 

1 

WeibuU 

1.46824 

1000.0 

100.0 

4795.784 

~ 

— 

1.0 

1.000 

3(4) 

WeibuU 

1.50438 

10.0 

4.17 

7.973 

4(5) 

WeibuU 

5.64760 

10.0 

1.39 

0.081 

5(6) 

WeibuU 

2.61249 

1000.0 

12.5 

26.442 

_ 6i7) _ 

WeibuU 

2.61249 

1000.0 

12.5 

26.442 

Setting  C. 

1 

WeibuU 

1.46824 

1000.0 

100.0 

4795.784 

~ 

~ 

1.0 

1.000 

■BBH 

WeibuU 

5.6476 

10.0 

1.39 

0.081 

4(5) 

WeibuU 

5.6476 

10.0 

1.39 

0.081 

WeibuU 

2.06810 

1000.0 

25.0 

160.798 

WeibuU 

2.61249 

1000.0 

12.5 

26.442 

Table  4.3.  Service  time  distribution  settings  for  the  closed  queueing  network 
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Design 

Model 

Transition  Matrix 

Service  Time  Setting 

Point 

Oi 

■EH 

Pi 

P2 

A. 

B. 

C. 

1 

X 

X 

X 

2 

X 

X 

X 

3 

X 

X 

X 

4 

X 

X 

X 

5 

X 

X 

X 

6 

X 

X 

X 

7 

X 

X 

X 

8 

X 

X 

X 

9 

X 

X 

X 

10 

X 

X 

X 

11 

X 

X 

X 

12 

X 

X 

X 

Table  4.4.  Design  points  for  closed  queueing  network 


4.3  Computer  Implementation 

All  simulation  models  are  programmed  in  the  simulation  language  SLAM  II  and 
implemented  on  a  Sun  SPARC  4  work  station.  Response  statistics,  internal  control 
variates,  and  external  control  variates  are  collected  and  calculated  using  FORTRAN  user- 
inserts  with  the  SLAM  11  programs.  For  the  open  queueing  network  of  section  3.3, 
analytical  control  variates  are  calculated  using  the  FORTRAN  insert.  The  Pascal  program 
ForQue  is  used  to  calculate  the  analytical  control  variates  for  the  closed  network. 
Controlled  responses,  estimated  variances,  and  confidence  intervals,  as  well  as  the 
statistical  comparisons  of  section  4.2.1  above,  are  calculated  using  the  MATLAB 
computer  language.  SLAM  11  and  FORTRAN  source  code  is  provided  in  appendix  A  for 
one  design  point  of  each  network.  This  also  includes  the  code  necessary  to  generate  the 
external  control  variates  for  that  design  point  (models  M3  and  O^).  Appendk  B  contains 
the  MATLAB  script  files  used  to  calculate  and  compare  controlled  response  performance. 
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Application  of  these  programs  at  a  design  point  is  a  three  step  process.  First  the 
model  under  study  (Mi,  Qi,  or  Qc)  is  replicated  1,000  times  using  the  SLAM  II  code. 
Following  the  end  of  each  replication,  the  FORTRAN  insert  collects  the  appropriate 
statistics  and  calculates  the  values  of  the  response  values,  internal  control  variates,  and  (in 
the  case  of  the  open  network)  analytical  control  variates.  These  values  are  then  put  into  a 
MATLAB  file  of  formatted  matrices  by  subroutines  in  the  FORTRAN  insert.  In  the  case 
of  the  closed  network,  the  appropriate  statistics  required  by  ForQue  to  calculate  the 
analytical  control  variates  are  put  in  a  text  file  by  the  FORTRAN  insert.  Following  all 
1,000  replications  the  text  file  is  read  by  ForQue  which  calculates  the  analytical  control 
variates  for  each  replication  and  puts  these  values  in  another  text  file.  Next,  the  external 
model  (Ms  or  Q3)  is  replicated  1,000  times  in  a  synchronized  manner.  Another 
FORTRAN  insert  collects  the  appropriate  statistics  following  each  replication  and  places 
the  external  control  variates  in  a  MATLAB  formatted  file.  Finally,  the  MATLAB  script 
file  reads  the  appropriate  MATLAB  and  text  files  and  calculates  uncontrolled  and 
controlled  responses  and  performs  the  necessary  comparisons. 

The  expected  values  of  the  steady-state  sojourn  time  for  the  open  network  and  the 
system  sojourn  time  and  CPU  utilization  fraction  for  the  closed  network  are  estimated  by 
making  25,000  replications  at  each  design  point.  While  these  values  are  estimates,  the 
associated  .90  confidence  intervals  are  sufficiently  tight  (less  than  0.5%  of  estimated  value 
in  all  cases)  to  make  good  estimates  of  coverage  and  Mean  Square  Error  for  comparison 
purposes. 
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4.4  Open  Queueing  Network  Results 

Analytical  control  variates  provide  significant  variance  and  confidence  interval 
width  reduction  on  the  estimates  of  customer  sojourn  time  at  each  of  the  3  design  points 
for  both  10  and  20  replications.  Since  confidence  interval  width  is  a  function  of  variance, 
and  of  primary  concern  in  most  simulation  studies,  only  confidence  interval  width 
reduction  is  presented.  Table  4.5  presents  confidence  interval  width  reduction  percentages 
achieved  for  the  analytical  and  external  control  variates  as  described  in  section  4.2. 
Results  for  the  “best”  combination  of  internal  control  variates  are  also  included.  The 
“best”  combination  is  defined  as  the  combination  of  internal  control  variates  that  produce 
the  most  confidence  interval  width  reduction. 

Achieved  coverage  percentages  for  uncontrolled  and  controlled  responses  are 
presented  in  table  4.6.  The  internal  control  variate  results  are  from  the  same  combination 
of  internal  control  variates  used  in  table  4.5.  Estimated  Mean  Square  Error  values  are 
listed  in  table  4.7. 


Traffic  Intensity 

Replications 

Confidence  Interval  Width  Reduction  (%)  | 

Analytical 
control  variates 

External  control 
variates 

Internal  control 
variates 

0.50 

10 

15.5 

29.7 

9.3 

0.50 

20 

19.7 

32.6 

13.9 

0.75 

10 

25.5 

34.4 

20.2 

0.75 

20 

28.5 

37.5 

24.1 

0.90 

10 

26.2 

38.7 

22.0 

0.90 

20 

28.7 

41.7 

25.1 

Table  4.5.  Confidence  interval  width  reduction  (%)  for  open  queueing  network 
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■ 

Replications 

Coverage  (%) 

Uncontrolled 

Analytical 

control 

variates 

External 

control 

variates 

Internal 

control 

variates 

0.50 

10 

86 

89 

88 

91 

0.50 

20 

86 

84 

90 

88 

0.75 

10 

87 

72 

84 

86 

0.75 

20 

88 

64 

80 

84 

0.90 

10 

87 

84 

84 

89 

0.90 

20 

92 

82 

82 

92 

Table  4.6.  Realized  coverage  (%)  for  open  queueing  network 


Traffic 

intensity 

Replications 

Mean  Square  Error 

Uncontrolled 

Anal5l:ical 

control 

variates 

External 

control 

variates 

Internal 

control 

variates 

0.50 

10 

0.222 

0.129 

0.093 

0.169 

0.50 

20 

0.132 

0.063 

0.037 

0.076 

0.75 

10 

16.003 

12.190 

6.551 

10.965 

0.75 

20 

7.641 

7.274 

3.141 

5.770 

0.90 

10 

360.694 

185.446 

126.482 

173.598 

0.90 

20 

186.145 

89.724 

75.283 

66.380 

Table  4.7.  Estimated  Mean  Square  Error  for  open  queueing  network 
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Several  interesting  results  become  evident  from  the  collected  data.  First,  in  terms  of 
confidence  interval  width  reduction,  external  control  variates  provide  the  best 
performance,  followed  by  analytical  and  then  internal  control  variates.  The  superior 
performance  of  the  external  control  variates  is  due  primarily  to  the  fact  that  the  customer 
arrival  sequence  is  exactly  the  same  for  Mi  and  M2.  Furthermore,  customers  then  proceed 
through  the  network  in  the  same  order  they  enter  the  network  and  experience  correlated 
service  times  in  only  three  service  centers  before  exiting  the  system.  Also,  the 
performance  of  all  of  the  control  variates  improves  as  trafBc  intensity  increases,  which 
corresponds  to  the  increasing  value  of  variance  of  the  uncontrolled  response  estimate. 

It  is  apparent  from  the  table  4.6  that  analytical  control  variates  may  be  inducing  bias 
into  the  controlled  estimate  of  customer  sojourn  time.  As  discussed  above  in  section  4.2, 

bias  can  be  introduced  in  the  process  of  estimating  .  Particularly  poor  coverage  is 
realized  for  the  traffic  intensity  of  0.75.  External  control  variates  also  indicate  a  possible 
problem  with  bias,  but  not  to  as  great  an  extent  as  that  of  the  analytical  control  variates. 
However,  internal  control  variate  realized  coverages  are  all  within  6  percent  of  the 
nominal. 

On  the  other  hand,  estimated  MSE  values  are  not  indicative  of  significant  bias  for 
analytical  control  variates.  In  all  cases  the  estimated  MSE  for  analytical  control  variates  is 
less  than  that  for  the  uncontrolled  sojourn  time  estimate.  It  should  be  noted  though,  that 
MSE  estimates  for  external  and  internal  control  variates  are  less  than  that  of  the  analytical 
control  variates,  except  for  one  case,  indicating  that  analytical  control  variates  are  not 
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performing  with  the  same  accuracy  as  that  of  the  other  types.  The  issue  of  bias  will  be 
explored  in  more  detail  for  the  closed  queuing  network  in  the  next  section. 

Another  interesting  observation  is  that  the  “best”  internal  control  variate  combination 
at  4  of  the  design  points  is  the  standardized  arrival  variable  and  at  the  other  2  design 
points  the  combination  of  the  standardized  arrival  variable  and  the  standardized  work 
variable  for  center  2  achieves  the  best  confidence  interval  width  reduction.  (It  should  be 
noted  that  the  2  design  points  where  2  internal  control  variates  achieve  the  most  variance 
reduction  consist  of  20  replications.  See  section  4.5.1  below  for  a  discussion  of  the  loss 
factor  as  a  probable  cause  of  this  behavior.)  As  a  point  of  reference,  table  4.8  lists  the 
confidence  interval  width  reduction  of  the  5  best  internal  control  variate  combinations  at 
the  design  point  of  traffic  intensity  equal  to  0.75  with  replication  size  of  20.  Due  to  the 
nature  of  the  network,  it  can  be  expected  that  customer  arrivals  should  drive  system 
performance  since  arrivals,  being  exponentially  distributed,  are  highly  variable.  The 
Weibull  service  times  are  not  as  variable.  The  performance  of  the  standardized  arrival 
variable,  and  the  implicit  high  correlation  between  it  and  customer  sojourn  time,  bears  out 
this  expectation.  For  the  open  queueing  network,  using  internal  control  variates  provides 
additional  insight  into  system  performance  that  analytical  and  external  control  variates 
don’t  necessarily  provide. 
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Internal  control 
variates 

Confidence 
interval  width 

Confidence 
interval  width 
reduction  (%) 

W2A 

6.71 

24.1 

A 

6.74 

23.7 

W^ 

6.92 

21.7 

W2W3A 

6.93 

21.6 

W1W2A 

7.02 

20.5 

Table  4.8.  Performance  of  5  best  internal  control  variate  combinations  for  open  queueing 
network  (traffic  intensity  =  0.75,  replications  =  20) 


4. 5  Closed  Queueing  Network  Results 

Results  for  the  closed  queueing  network  are  presented  below.  A  discussion  of 
confidence  interval  width  reduction  is  presented  followed  by  a  discussion  on  estimation  of 
bias. 

4.5.1  Confidence  Interval  Width  Reduction.  Analytical  control  variates  provide 
significant  confidence  interval  reduction  on  estimates  for  both  system  sojourn  time  and 
CPU  utilization  for  the  closed  queueing  network.  Across  the  range  of  all  design  points, 
except  for  a  few  isolated  cases,  analytical  control  variate  performance  is  similar  to  that  of 
external  and  internal  control  variates.  Confidence  interval  vridth  reductions,  as  a 
percentage  of  the  uncontrolled  estimated  confidence  interval  are  provided  in  tables  4.9 
through  4.12.  System  sojourn  time  results  are  listed  in  tables  4.9  and  4.10  while  results 
for  CPU  utilization  are  enumerated  in  tables  4.11  and  4.12.  Results  are  included  for 
analytical  control  variates,  external  control  variates,  and  internal  control  variates.  The 
internal  control  variates  form  the  combination  of  standardized  work  variables  and 
standardized  routing  variables,  that  produces  the  greatest  reduction  in  confidence  interval 
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width.  It  should  be  noted  that  for  one  design  point  (when  measuring  CPU  utilization  for 
model  Oc)  external  control  variates  fail  to  produce  any  variance  reduction.  Where  this 
occurs,  the  symbol  is  used  in  the  table. 

Of  particular  interest  is  the  performance  of  external  and  analytical  control  variates  for 
model  Oc-  External  control  variates  fail  to  provide  the  same  level  of  confidence  interval 
width  reduction,  particularly  for  CPU  utilization,  for  model  Qc  as  they  do  for  model  Qi. 
Recall  that  in  model  Qc,  only  5  customers  at  a  time  are  allocated  a  portion  of  main 
memory,  and  therefore  access  to  the  CPU  and  mass  storage  devices  (centers  3  through  7). 
On  the  other  hand,  model  O3  has  an  unconstrained  subnetwork  and  an  unlimited  number 
of  customers  may  access  the  CPU  and  storage  devices.  Certainly  the  variance  in  external 
control  variate  performance  is  due  to  this  dissimilarity  in  model  structure.  This  difference 
causes  the  sequence  of  customers  visiting  a  particular  service  center  to  differ  even  more 
than  when  model  Oi  is  the  model  under  study.  Under  these  conditions,  the  common 
random  numbers  lose  even  more  synchronization  and  the  system  responses  from  each 
model  are  not  as  highly  correlated. 

Analytical  control  variates,  on  the  other  hand,  continue  to  perform  at  the  same  level 
for  Qc-  Although  the  underlying  model,  O2,  is  also  unconstrained,  the  model  relies  only 
on  the  mean  responses  of  system  parameters  from  Qc-  Therefore,  the  same  conditions 
that  cause  longer  sojourn  times  or  greater  CPU  utilization  levels  in  model  Qc  will  do  the 
same  in  Qx-  Hence,  correlation  from  replication  to  replication  is  maintained  and  analytical 
control  variates  continue  to  reduce  variance  for  Qc. 
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Model  0i  -  System  sojourn  time 


Service  time 
setting 

Transition 

probability 

matrix 

Confidence  interval  width  reduction  (%) 
(Replications  =  10(20)) 

Analytical 
control  variate 

External  control 
variate 

Internal  control 
variates 

A. 

Pi 

51.8(52.5) 

44.9(56.9) 

Pi 

42.5(43.7) 

B. 

Pi 

46.1(46.9) 

63.0(64.2) 

56.8(66.4) 

Pi 

47.1(48.4) 

42.0(52.4) 

C. 

Pi 

52.6(55.7) 

Pi 

18.8(21.9) 

47.6(50.8) 

57.6(64.3) 

Table  4.9.  Confidence  interval  width  reduction  for  system  sojourn  time  for  model  Qi 


Model  Qc  -  System  sojourn  time 

Service  time 
setting 

Transition 

probability 

matrix 

Confidence  interval  Avidth  reduction  (%) 
(Replications  =  10(20)) 

Analytical 
control  variate 

External  control 
variate 

Internal  control 
variates 

A. 

Pi 

WW-fM 

IKXESSSHi 

Pi 

28.5(32.3) 

32.2(46.0) 

B. 

Pi 

■cVXMlM 

59.8(68.0) 

Pi 

56.8(57.2) 

20.4(22.4) 

40.4(46.8) 

C. 

Pi 

19.9(24.4) 

37.9(49.7) 

Pi 

55.8(57.9) 

7.4(10.6) 

63.3(68.6) 

Table  4. 10.  Confidence  interval  width  reduction  for  system  sojourn  time  for  model  Oc 
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1  Moc 

el  Oi  -  CPU  utilization 

Service  time 
setting 

Transition 

probability 

matrix 

Confidence  interval  width  reduction  (%) 
(Replications  =  10(20)) 

Analytical 
control  variate 

External  control 
variate 

Internal  control 
variates 

A. 

Pi 

54.2(55.6) 

45.0(46.8) 

39.0(47.1) 

Pi 

57.4(60.5) 

31.7(54.3) 

B. 

Pi 

6S.7(70.6) 

61.8(63.3) 

77.8(81.1) 

Pi 

54.4(57.4) 

49.8(52.0) 

32.2(47.8) 

C. 

Pi 

39.7(42.4) 

54.0(55.8) 

18.8(34.8) 

Pi 

41.6(44.9) 

49.0(51.6) 

80.1(83.6) 

Table  4. 1 1 .  Confidence  interval  width  reduction  for  CPU  utilization  for  model  Qi 


Mod 

el  Qc  -  CPU  utilization 

Service  time 
setting 

Transition 

probability 

matrix 

Confidence  interval  width  reduction  (%) 
(Replications  =  10(20)) 

Analytical 
control  variate 

External  control 
variate 

Internal  control 
variates 

A. 

Pi 

37.2(38.7) 

12.5(16.2) 

16.6(34.9) 

Pi 

41.2(43.3) 

5.6(11.0) 

9.4(30.9) 

B. 

Pi 

63.8(65.4) 

59.8(65.5) 

Pi 

29.4(32.5) 

~(3.6) 

16.8(32.4) 

C. 

Pi 

42.2(44.3) 

4.9(10.5) 

30.1(42.3) 

Pi 

68.3(69.3) 

5.3(9.7) 

67.3(71.3) 

Table  4.12.  Confidence  interval  width  reduction  for  CPU  utilization  for  model  Qc 
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Replication  length  seems  to  influence  internal  control  variate  performance  more  than 
that  of  the  other  control  variates.  Analytical  and  external  control  variates  produce  only 
slightly  more  confidence  interval  width  reduction  when  20  replications  are  used  versus  10 
replications.  However,  for  internal  control  variates,  the  increased  number  of  replications 
produce  significantly  more  confidence  interval  width  reduction  than  when  only  10 
replications  are  used.  On  average,  internal  control  variates  produce  11.2  percent  more 
confidence  interval  width  reduction  for  20  versus  10  replications  across  all  design  points. 
Analytical  and  external  control  variates  only  produce  1.9  percent  and  2.7  percent  more 
reduction  for  20  replications. 

The  poorer  performance  of  the  internal  control  variates  at  10  replications  can  be 
explained  by  the  loss  factor  discussed  in  section  1.2.1.  Recall  that  the  variance  of  a 
controlled  point  estimate  is  given  by 

Far[y(y0)]  =  (l  -  )Far[F]  (4.6) 

where  n  is  the  number  of  replications,  q  is  the  number  of  control  variates,  and  is  the 
coefficient  of  correlation  between  the  random  variable  Y  and  the  vector  of  control  variates, 
and  the  loss  factor  is  («  -  2)/(n  -  ^  -  2) .  For  the  closed  queueing  network,  the  average 

number  of  internal  control  variates  in  the  “best”  combination  is  3.21  for  10  replications 
and  4.46  for  20  replications.  Calculating  the  loss  factor  for  3  control  variates  yields  1.6 
for  10  replications  and  1.2  for  20  replications.  For  4  control  variates  the  loss  factor  is  2 
and  1.29  respectively.  The  loss  factors  for  both  analytical  and  external  control  variates 
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remain  fixed  at  1.14  and  1.05  for  10  and  20  replications,  since  each  is  a  single  control 
variate.  Obviously,  at  10  replications  the  loss  factor  has  a  higher  detrimental  effect  on 
variance  reduction  for  the  larger  combinations  of  internal  control  variates  than  it  does  on 
the  analytical  and  external  control  variate. 

As  with  the  open  queueing  network,  internal  control  variates  provide  additional 
insight  into  system  performance  at  each  design  point.  By  studying  the  reduction  achieved 
by  the  different  combinations  of  internal  control  variates,  critical  system  characteristics 
become  apparent.  For  example,  table  4.13  lists  the  confidence  interval  width  reduction 
achieved  by  the  5  best  combinations  of  internal  control  variates  for  CPU  utilization.  The 
design  point  is  model  Qi,  with  service  time  setting  C.,  probability  transition  matrix  P2,  and 
replications  equal  to  20.  Obviously,  the  proportion  of  customers  routed  to  center  5,  and 
the  deviation  from  the  mean  service  time  at  center  5  play  a  significant  role  in  system 
performance. 


Internal  control 

Confidence  interval 

variate  combination 

width  reduction  (%) 

W2WSR5 

83.6 

W1W2W5R5 

83.4 

83.2 

W2W^xR, 

83.2 

W2WJilR^R, 

82.6 

Table  4.13.  Confidence  interval  width  reduction  (%)  for  CPU  utilization.  Design  point: 
Qi,  service  time  setting  C.,  F2,  and  20  replications. 


The  results  of  analytically  controlled  responses  themselves  cannot  provide  the  same 
insight.  However,  the  mean  value  analysis  algorithm  used  by  ForQue,  the  program  used 
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to  calculate  the  analytical  control  variate,  provides  more  information  about  the  network 
than  just  system  sojourn  time  and  CPU  utilization.  Given  the  good  performance  of  the 
analytical  control  variates  for  both  performance  measures  across  all  design  points,  the 
output  of  ForQue  must  be  highly  correlated  to  the  output  of  both  Oi  and  Qc-  In  a  sense, 
the  variance  reduction  performance  of  the  analytical  control  variates  provides  validation 
for  the  analytical  model’s  ability  to  predict  model  Oi  and  Oc  behavior.  For  example,  table 
4.14  is  the  output  from  ForQue  using  the  mean  values  of  the  system  parameters  of  the 
same  design  point  described  above  for  table  4.13.  The  extremely  high  utilization  fraction 
and  significantly  longer  response  time  at  center  5  predicted  by  ForQue  point  to  the  same 
inferences  about  system  performance  as  those  made  in  the  previous  section  by  analyzing 
the  internal  control  variates. 


***  FORQUE  PERFORMANCE  REPORT  *** 
Number  of  Customers  =  25 


sta¬ 

Nmbr 

Average 

Visit 

Through 

Queue 

Respons 

Utiliz 

tion 

Chls 

Svc  Tm 

Ratio 

-put 

Length 

Time 

-ation 

1 

25 

100.0000 

1.0000 

0.1317 

13.1747 

100.0000 

13.1747 

2 

1 

1.0000 

4.0000 

0.5270 

1.0940 

2.0759 

0.5270 

3 

1 

1.3900 

1.2000 

0.1581 

0.2810 

1.7775 

0.2198 

4 

1 

1.3900 

1.2000 

0.1581 

0.2810 

1.7775 

0.2198 

5 

1 

25.0000 

0.3000 

0.0395 

9.2069 

232.9453 

0.9881 

6 

1 

12.5000 

0.3000 

0.0395 

0.9624 

24.3501 

0.4940 

Sojourn  time  =  89.7583 


Table  4. 14.  Output  from  ForQue  computer  program.  Network  settings:  Q2,  service  time 
setting  C.,  transition  probability  matrix  P2. 
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4.5.2  Estimation  of  Bias.  Realized  confidence  interval  coverage  for  uncontrolled 
and  controlled  response  estimates  are  enumerated  in  tables  4.15  through  4.18.  Nominal 
coverage  is  90  percent.  Due  to  the  similarity  in  coverage  achieved,  only  results  for  20 
replications  are  provided.  Internal  control  variate  coverage  corresponds  to  the  same 
internal  control  variate  combinations  reported  above  for  confidence  interval  width 
reduction.  As  with  the  open  queuing  network,  analytical  control  variates  seem  to  induce 
bias  into  the  controlled  estimates  of  system  sojourn  time  and  CPU  utilization.  Internal  and 
external  control  variate  coverage  results  indicate  little,  if  any,  problem  with  bias  for  their 
controlled  estimates. 

Estimated  Mean  Square  Error  values  are  listed  in  tables  4.19  through  4.22.  Despite 
the  indications  of  bias  from  the  realized  coverage  percentages,  estimated  MSE  for  the 
analytical  control  variates  is  substantially  less  than  that  of  the  uncontrolled  responses  in 

almost  every  case.  This  indicates  that  both  and  are  more 

accurate,  on  the  average,  in  estimating  system  sojourn  time  and  CPU  utilization 
thanr(20)  andf/(20).  However,  some  combination  of  the  reduced  confidence  interval 
width  and  real  bias  is  preventing  analytical  control  variates  from  achieving  reasonable 
levels  of  realized  coverage.  Other  indicators  need  to  be  examined  to  determine  what  is 
affecting  the  coverage  values. 
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Model  0\  -  System  sojourn  time 


Service  time 
setting 

Transition 

probability 

matrix 

Coverage  (%) 

Uncontrolled 

response 

Analytical 

control 

variate 

External 

control 

variate 

Internal 

control 

variate 

A. 

Pi 

92 

68 

90 

86 

Pi 

86 

64 

90 

84 

B. 

Pi 

86 

78 

94 

88 

Pi 

84 

62 

92 

88 

C. 

Pi 

94 

72 

82 

88 

Pi 

88 

88 

88 

90 

Table  4. 15.  Realized  coverage  for  system  sojourn  time  for  model  Q\  (replications  =  20) 


Model  Oc  -  System  sojourn  time 


Service  time 
setting 

Transition 

probability 

matrix 

Coverage  (%) 

Uncontrolled 

response 

Analytical 

control 

variate 

External 

control 

variate 

Internal 

control 

variate 

A. 

Pi 

90 

66 

92 

80 

Pi 

88 

50 

94 

90 

B. 

Pi 

90 

78 

90 

94 

Pi 

94 

38 

96 

88 

C. 

Pi 

98 

1  52 

96 

96 

Pi 

88 

1  92 

86 

88 

Table  4. 16.  Realized  coverage  for  system  sojourn  time  for  model  Qc  (replications  =  20) 
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Model  Oi  -  CPU  utilization 


Service  time 
setting 

Transition 

probability 

matrix 

Coverage  (%) 

Uncontrolled 

response 

Analytical 

control 

variate 

External 

control 

variate 

Internal 

control 

variate 

A. 

Pi 

90 

68 

92 

88 

Pi 

90 

62 

84 

88 

B. 

Pi 

88 

82 

80 

92 

Pi 

92 

58 

92 

90 

C. 

Pi 

92 

54 

82 

88 

Pi 

86 

92 

80 

86 

Table  4.17.  Realized  coverage  for  CPU  utilization  for  model  Oi  (replications  =  20) 


Model  Qc  -  CPU  utilization 

Service  time 
setting 

Transition 

probability 

matrix 

Coverage  (%) 

Uncontrolled 

response 

Analytical 

control 

variate 

External 

control 

variate 

Internal 

control 

variate 

A. 

Pi 

90 

88 

94 

86 

Pi 

84 

82 

84 

86 

B. 

Pi 

92 

92 

94 

88 

Pi 

88 

90 

88 

84 

C 

Pi 

90 

68 

88 

96 

Pi 

90 

78 

86 

92 

Table  4.18.  Realized  coverage  for  CPU  utilization  for  model  Qc  (replications  =  20) 
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Model  Oi  -  System  sojourn  time 


Service  time 
setting 

Transition 

probability 

matrix 

Estimated  Mean  Square  Error 

Uncontrolled 

response 

Analytical 

control 

variate 

External 

control 

variate 

Internal 

control 

variate 

A. 

Pi 

3.146 

1.609 

0.713 

0.629 

Pi 

2.976 

1.635 

0.511 

0.842 

B. 

Pi 

22.530 

6.509 

1.993 

2.070 

Pi 

5.244 

2.602 

0.590 

1.070 

C. 

Pi 

5.671 

5.881 

1.473 

2.480 

Pi 

38.641 

18.654 

8.239 

3.938 

Table  4.19.  Estimated  Mean  Square  Error  for  system  sojourn  time  for  model  0i 
(replications  =  20) 


Model  Qc  -  System  sojourn  time 


Service  time 
setting 

Transition 

probability 

matrix 

Estimated  Mean  Square  Error 

Uncontrolled 

response 

Analytical 

control 

variate 

External 

control 

variate 

Internal 

control 

variate 

A. 

Pi 

5.830 

2.635 

1.390 

1.623 

Pi 

6.387 

3.876 

1.730 

1.413 

B. 

Pi 

25.902 

5.902 

7.755 

1.679 

Pi 

8.584 

8.006 

4.579 

2.356 

C. 

Pi 

15.605 

16.121 

8.027 

2.623 

Pi 

54.500 

7.247 

35.259 

4.264 

Table  4.20.  Estimated  Mean  Square  Error  for  system  sojourn  time  for  model  Qc 
(replications  =  20) 
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Model  Qi  -  CPU  utilization 


Service  time 
setting 

Transition 

probability 

matrix 

Estimated  Mean  Sc 

luare  Error  (*  10'^) 

Uncontrolled 

response 

Analytical 

control 

variate 

External 

control 

variate 

Internal 

control 

variate 

A. 

Pi 

6.44 

3.55 

2.36 

2.09 

Pi 

9.08 

4.19 

2.88 

2.59 

B. 

Pi 

10.72 

1.08 

1.77 

0.48 

Pi 

5.82 

3.10 

1.26 

1.49 

C. 

Pi 

18.10 

17.96 

5.88 

7.33 

Pi 

26.70 

6.46 

8.84 

0.55 

Table  4.21.  Estimated  Mean  Square  Error  (*  10'^)  for  CPU  utilization  for  model  Qi 
(replications  =  20) 


Model  Qc  -  CPU  utilization 

Service  time 
setting 

Transition 

probability 

matrix 

Estimated  Mean  Sc 

luare  Error  (*  10'^) 

Uncontrolled 

response 

Analytical 

control 

variate 

External 

control 

variate 

Internal 

control 

variate 

A. 

Pi 

4.82 

2.61 

3.18 

2.72 

Pi 

11.10 

3.86 

9.02 

3.84 

B. 

Pi 

7.35 

0.66 

4.11 

0.70 

Pi 

5.25 

2.15 

4.81 

2.37 

C 

Pi 

18.10 

12.90 

14.70 

3.66 

Pi 

18.60 

3.11 

13.20 

1.77 

Table  4.22.  Estimated  Mean  Square  Error  (*  10'*)  for  CPU  utilization  for  model  Qc 
(replications  =  20) 
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Another  indication  of  the  coverage  problem  for  analytical  control  variates  is  presented 
in  tables  4.23  and  4.24.  These  tables  list  the  estimated  mean  for  each  design  point  found 
using  25,000  replications  and  the  average  controlled  response  estimated  using  analytical 

control  variates  and  )•  For  system  sojourn  time,  the  analytically 

controlled  response  is  lower  than  the  estimated  response  at  every  design  point.  The 
situation  is  reversed  for  CPU  utilization  where  the  analytically  controlled  estimate  is 
consistently  higher,  with  the  exception  of  2  design  points.  Both  of  these  points 
correspond  to  the  same  service  time  setting  (C.)  and  transition  probability  matrix  {Pi). 

Tables  4.25  and  4.26  compare  the  realized  estimates  for  system  sojourn  time  and 
CPU  utilization  with  the  mean  values  obtained  at  each  design  point  using  ForQue.  Note 
that  the  mean  value  analysis  of  the  exponential  network  provided  by  ForQue  is 
consistently  higher  for  system  sojourn  time  and  consistently  lower  for  that  of  CPU 
utilization  than  the  estimates  obtained  with  25,000  replications  for  Qx.  These  relationships 
are  directly  opposite  of  that  for  the  analytically  controlled  response  and  the  estimated 
means.  These  relationships  could  have  some  bearing  on  the  problem. 
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Table  4.23.  Estimated  system  sojourn  time  for  25,000  replication  vs  analytical  controlled 
response  for  20  replications 


CPU  utilization 


Transition 

probability 

matrix 


C(25,000) 

Percentage 

difference 

0.903 

0.907 

0.44 

0.755 

0.761 

0.79 

0.663 

0.664 

0.15 

0.686 

0.691 

0.73 

0.850 

0.861 

1.29 

0.539 

0.536 

-0.56 

0.866 

0.868 

0.23 

0.708 

0.712 

0.56 

0.618 

0.618 

0.00 

0.612 

0.615 

0.49 

0.775 

0.784 

3.74 

0.516 

0.512 

-0.77 

Table  4.24.  Estimated  CPU  utilization  for  25,000  replication  vs  analytical  controlled 
response  for  20  replications 
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System  sojourn  time 

Model 

Service  time 
setting 

Transition 

probability 

matrix 

7’(25,000) 

Mean  value 
analysis 
(ForQue) 

01 

A. 

Pi 

38.9 

41.2 

P2 

32.9 

37.2 

B. 

Pi 

88.5 

91.1 

P2 

46.3 

51.0 

C. 

Pi 

47.8 

52.3 

P2 

87.7 

89.8 

Table  4.25.  Comparison  of  estimated  system  sojourn  time  and  mean  value  analysis 


CPU  utilization 

Model 

Transition 

probability 

matrix 

C/(25,000) 

Mean  value 
analysis 
(ForQue) 

01 

A 

Pi 

0.903 

0.885 

P2 

0.755 

0.729 

B. 

Pi 

0.663 

0.654 

P2 

0.686 

0.662 

C. 

Pi 

0.850 

0.821 

P2 

0.539 

0.527 

Table  4.26.  Comparison  of  estimated  CPU  utilization  and  mean  value  analysis 


To  get  a  better  idea  of  exactly  how  this  induced  bias  is  affecting  the  analytically 
controlled  estimates,  the  following  figures  are  provided.  Each  figure  depicts  the 
estimated  responses  and  associated  confidence  intervals  for  all  experiments  at  a  single 
design  point.  For  system  sojourn  time,  model  Oi  with  service  time  setting  B.,  transition 
probability  matrix  P2,  and  20  replications  (50  experiments)  is  presented.  The  50  diamond 
symbols  in  each  figure  represent  the  50  point  estimates  made  at  this  design  point  using  20 
replications.  The  bracketed  lines  above  and  below  each  diamond  represent  the  width  of 
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the  estimated  confidence  interval  for  the  associated  20  replication  design  point.  For 
reference  the  estimated  mean  found  using  25,000  replications  is  represented  by  the 
horizontal  line  across  each  figure.  Figure  4. 1  presents  the  values  obtained  by  7(20)  and 
figures  4.2  and  4.3  compare  and  T{j3)^  for  the  same  experiments.  Figure 

4.2  clearly  shows  how  consistently  underestimates  system  sojourn  time  while 

appears  to  be  evenly  scattered  about  the  estimated  mean.  Although  the 
analytically  controlled  estimates  are  mostly  lower  than  the  mean,  on  average  they  are 
closer  than  the  uncontrolled  estimate,  hence  the  lower  estimated  MSE.  T(^)^  on  the 
other  hand,  appears  to  exhibit  no  evidence  of  bias. 


0  r(20) 

I  Confidence  interval 
: - Mean  =  46.3 


Figure  4. 1  7(20)  experimental  results  with  estimated  confidence  interval.  Design  point: 
Qi,  service  time  setting  B.,  transition  probability  matrix  P2.  (Each  point 
represents  one  point  estimate  and  confidence  interval  realization.) 
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0 

I  Confidence  interval 
- Mean  =  46.3 

Figure  4.2.  experimentaJ  results  with  estimated  confidence  interval.  Design 

point;  Qi,  service  time  setting  B.,  transition  probability  matrix  P2.  (Each 
point  represents  one  point  estimate  and  confidence  interval  realization.) 


I  Confidence  interval 
- Mean  =  46.3 


Figure  4.3.  T(P)^  experimental  results  with  estimated  confidence  interval.  Design 
point;  Oi,  service  time  setting  B.,  transition  probability  matrix  P2.  (Each 
point  represents  one  point  estimate  and  confidence  interval  realization.) 
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Finally,  model  Qc  with  service  time  setting  C.,  probability  transition  matrix  P2,  and  20 
replications  is  used  to  create  the  same  figures  for  CPU  utilization.  Figure  4.4  shows  the 
experimental  results  for  t/(20),  and  figures  4.5  and  4.6  present  the  experimental  results 

for  and  .  The  internal  control  variates  used  are  the  same  as  reported 

above  for  confidence  interval  width  reduction.  As  above,  the  analytically  controlled 
response  is  consistently  on  one  side  of  the  estimated  mean  while  consistently  closer  to  the 
mean  than  the  uncontrolled  response.  Similarly,  the  internal  control  variates  are  evenly 
scattered  about  the  mean  and  exhibit  no  evidence  of  bias. 

The  above  figures  coupled  with  the  data  on  realized  coverage  and  the  consistent 
difference  between  the  estimated  mean  and  analytically  controlled  responses  point  to  a  real 
problem  with  bias  for  the  analytical  control  variates  in  this  study.  The  mechanism(s) 
causing  the  bias  in  the  analytically  controlled  response  estimates  is(are)  not  totally  clear. 
Means  of  possibly  reducing  the  bias  are  discussed  in  the  following  chapter. 


0  C/(20) 

I  Confidence  interval 
- Mean  =  .5158 


Figure  4.4.  U (20)  experimental  results  with  estimated  confidence  interval.  Design  point: 
Qc,  service  time  setting  C.,  transition  probability  matrix  P2.  (Each  point 
represents  one  point  estimate  and  confidence  interval  realization.) 
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I  Confidence  interval 
- Mean  =  .5158 


Figure  4.5.  experimental  results  with  estimated  confidence  interval.  Design 

point;  2c,  service  time  setting  C.,  transition  probability  matrix  P2.  (Each 
point  represents  one  point  estimate  and  confidence  interval  realization.) 


I  Confidence  interval 
- Mean  =  .5158 


Figure  4.6.  U experimental  results  with  estimated  confidence  interval.  Design 
point;  2c,  service  time  setting  C.,  transition  probability  matrix  Pj.  (Each 
point  represents  one  point  estimate  and  confidence  interval  realization.) 
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5.  Conclusions  and  Recommendations  for  Future  Research 


5.1  Introduction 

The  motivation  for  this  research  was  to  determine  if  a  hybrid  type  of  control  variate, 
called  an  analytical  control  variate  in  this  study,  could  effectively  reduce  the  width  of  point 
estimate  confidence  intervals  of  replicative  simulation  studies  while  avoiding  some  of  the 
limitations  of  internal  and  external  control  variates.  In  terms  of  confidence  interval  width 
reduction,  the  experimental  results  indicate  that  analytical  control  variates  are  quite 
successful  for  the  networks  studied.  Unfortunately  point  estimates  found  using  anal5l;ical 
control  variates  appear  to  have  significant  levels  of  bias.  However,  there  may  be 
occasions  where  the  bias  is  not  a  problem.  For  example,  analytical  control  variates  could 
be  useful  when  assessing  the  percentage  change  in  performance  when  system  parameters 
are  changed.  The  following  sections  contain  observations  on  the  differences  in  application 
of  anal5^ical  control  variates  versus  external  and  internal  control  variates.  Additionally,  a 
discussion  on  possible  remedies  to  the  bias  problem  is  presented.  Recommendations  for 
future  research  end  this  chapter. 

5.2  Study  Time  Efficiencies 

Analytical  control  variates  outperform  external  control  variates  in  this  study  in  terms 
of  time  efficiency.  For  example,  to  generate  the  external  control  variates  for  the  1,000 
replications  at  one  design  point,  1,000  replications  of  the  external  simulation  model  must 
be  made  as  well.  The  similarity  of  the  two  models  translates  into  approximately  equal 
computer  run  times  for  each  model.  Of  course,  in  some  computers,  such  as  the  Sun 
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SPARC  4  work  stations  used  in  this  study,  these  replications  can  be  run  at  the  same  time. 
However  parallel  processing  of  this  type  increases  average  run  time  for  each  model. 
Obviously,  only  truly  significant  gains  in  variance  reduction  over  other  control  variates 
would  make  external  control  variates  worth  while.  If  there  is  enough  time  to  make  the 
additional  external  simulations,  there  is  time  to  make  that  many  more  replications  of  the 
original  model  and  apply  internal  or  analytical  control  variates  to  the  additional 
replications. 

Computer  run  times  for  analytical  control  variates  are  virtually  instantaneous  in 
comparison.  For  the  open  system  studied,  analytical  control  variates  are  calculated  in  the 
FORTRAN  routine  used  to  collect  replication  statistics,  so  analytical  control  variates 
require  only  a  few  simple  calculations  following  each  replication.  Approximate  run  times 
for  the  closed  queueing  system  studied  are  15  minutes  for  1,000  replications  on  a  Sun 
SPARC  4  work  station.  The  ForQue  computer  program  can  generate  the  1,000  analytical 
control  variates  for  these  replications  in  less  than  10  seconds  on  the  same  computer 
system. 

The  generation  of  internal  control  variates  is  virtually  instantaneous  as  well.  In  both 
models  studied,  the  standardized  arrival,  work,  and  routing  variables  are  calculated  in  the 
same  FORTRAN  routine  that  collects  the  statistics  necessary  to  compute  them.  Although 
the  generation  of  the  internal  control  variates  is  as  rapid  as  that  for  analytical  control 
variates,  application  of  internal  control  variates  may  not  be  as  quick. 

Unlike  analytical  (and  external)  control  variates,  some  effort  (and  time)  must  be 
spent  in  choosing  the  “best”  subset  of  all  possible  internal  control  variates.  First  of  all,  not 
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all  of  the  internal  control  variates  generated  are  sufficiently  correlated  with  the  system 
output  to  produce  any  variance  reduction.  Secondly,  as  pointed  out  in  section  1.2.1,  at 
some  point  adding  additional  control  variates  that  are  correlated  with  the  system  output 
cannot  overcome  the  effect  of  the  loss  factor  in  the  estimated  variance.  Research  has  been 
conducted  on  the  best  way  to  select  an  optimum,  or  close  to  optimum,  set  of  internal 
control  variates.  Two  approaches  are  presented  in  Nelson  (1987)  (see  section  2.2.2). 
One  approach  involves  step-wise  regression  and  the  other  proposes  the  computation  of  a 
marginal  improvement  ratio.  Bauer  and  Wilson  (1992)  develop  a  method  for 
multiresponse  simulations  that  minimizes  the  mean-square  confidence-region  volume. 

One  obvious  time  advantage  to  internal  control  variates  over  both  analytical  and 
external  control  variates  is  that  no  analytical  model  that  approximates  the  behavior  of  the 
system  under  study  is  necessary.  Depending  on  the  system  simulated  finding  or  creating 
an  appropriate  analytical  model  can  be  a  daunting  problem,  however  there  are  many  useful 
analytical  models  in  existence  and  research  continues  in  the  development  of  more.  The 
mean  value  analysis  algorithm  used  for  the  closed  network  in  this  study  is  just  one 
example. 

From  this  study,  it  appears  that  if  an  adequate  analytical  model  is  available,  analytical 
control  variates  provide  an  advantage  in  terms  of  study  completion  times.  Further,  if  a 
particular  simulation  model  is  used  regularly  to  conduct  many  studies,  time  spent 
developing  or  finding  an  appropriate  analytical  model  could  pay  time  dividends  in  the 
future.  Otherwise,  internal  control  variates  are  the  only  option  available. 
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5.3  Effectiveness  over  Different  Model  Configurations 

As  the  synchronization  of  common  random  numbers  between  the  model  under  study 
and  the  external  model  decreased,  the  effectiveness  of  the  external  control  variates 
decreased.  Obviously,  the  correlation  between  the  outputs  of  the  two  models  depends  on 
the  correlation  of  the  input  random  variables.  As  models  increase  in  complexity, 
synchronicity  becomes  problematic.  Furthermore,  the  service  time  distributions  for  the 
models  under  study  in  this  research  effort  (Weibull)  have  the  characteristic  of  being 
generated  using  the  inverse  transform  method.  This  produces  a  one  to  one  relationship 
between  the  Weibull  random  variable  and  the  uniform  random  number  stream  used  to 
generate  them.  This  characteristic  made  the  level  of  synchronization  achieved  in  this  study 
possible.  Many  other  random  variable  distributions  that  may  be  required  to  adequately 
model  some  other  situation  can’t  be  generated  using  the  inverse  transform  method  and 
therefore  require  several  random  numbers  to  generate  one  random  variable  [13J. 
Achieving  any  random  number  synchronization  under  those  circumstances  becomes 
extremely  difficult  and  can  increase  computer  run  time  as  well. 

Analytical  control  variates  rely  only  on  the  correlation  between  one  replication  of  the 
model  under  study  and  the  result  of  an  analytical  model  calculation  using  the  observed 
mean  values  of  the  input  random  variables  for  that  replication.  As  long  as  the  analytical 
model  responds  in  a  correlated  manner,  variance  reduction  occurs.  In  this  research  effort, 
analytical  control  variates  perform  essentially  the  same  for  all  models  and  design  points 
considered  irregardless  of  their  departure  from  the  analytical  model.  Both  Tomick  (1988) 
and  Dietz  and  Harmonosky  (1989)  reported  similar  results. 
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Internal  control  variates  don’t  depend  on  some  external  approximation  to  the  model 
under  study  and  model  configuration  had  little  effect  on  their  performance.  Traffic 
intensities  and  utilization  rates  seemed  to  have  the  most  effect  on  internal  control  variate 
performance.  Variance  reduction  increases,  in  this  study,  for  internal  control  variates 
when  the  model  under  study  has  one  or  more  service  centers  with  high  traffic  intensity  or 
utilization  rates  inducing  high  correlation  between  the  internal  control  variates  that 
correspond  to  that  service  center  and  the  system  output. 

5.4  System  Insight 

As  discussed  previously  in  section  4.5,  the  variance  reduction  achieved  for  the 
different  internal  control  variates  provides  insight  into  the  behavior  of  the  network  under 
study.  The  internal  control  variates  that  produce  the  most  variance  reduction  indicate  that 
behavior  (service  times,  routing  probabilities)  at  the  service  centers  associated  with  them 
has  the  highest  correlation  to  system  performance. 

External  and  analytical  control  variates  can’t  provide  this  same  post-simulation 
analysis.  However,  the  analytical  models  used  to  produce  the  analytical  control  variates 
can  predict  similar  results  if  the  model  provides  system  parameter  values  (see  section  4.5). 
Obviously  these  predictions  are  only  as  good  as  the  analytical  model  approximates  the 
system  under  study  behavior.  As  discussed  in  section  4.5,  the  variance  reduction 
produced  by  a  particular  analytical  control  variate  is  a  function  of  the  correlation  between 
the  analytical  model  and  the  system  under  study.  Significant  variance  reduction  across  a 
range  of  system  parameters  indicates  that  the  analytical  model  output  tracks  with  the 
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output  of  the  model  under  study.  Given  this  variance  reduction  validation,  the  analytical 
model  could  be  used  by  analysts  to  quickly  provide  decision  makers  with  useful  answers 
when  there  is  not  enough  time  to  make  simulation  studies.  The  analytical  models  could  be 
used  in  much  the  same  way  response  surface  meta-models  are  now  used. 

Analytical  control  variates  can  also  provide  pre-simulation  analysis  that  neither  of  the 
other  two  types  of  control  variates  can.  The  analytical  model  used  to  generate  analytical 
control  variates  can  also  be  used  to  find  excellent  starting  points  for  studies  when  different 
configurations  or  resources  are  being  compared.  Neither  external  nor  internal  control 
variates  can  provide  the  analyst  with  the  same  information  before  any  replications  are 
made.  Given  a  good  starting  point,  study  times  for  large  projects  can  be  significantly 
reduced. 

5.5  Loss  Factor  and  Confidence  Interval  Width 

Discussion  in  Chapter  1  pointed  out  that  the  variance  estimator  when  control  variates 
are  used  includes  a  loss  factor  of  («-2)/(«-^-2)  where  n  is  the  number  of  replications  and  q 
is  number  of  control  variates.  Further,  the  degrees  of  freedom  for  the  associated  ^-statistic 
used  to  estimate  a  confidence  interval  is  equal  to  n-q-\.  Of  course,  if  an  additional  control 
variate  is  sufficiently  correlated  to  the  response  statistic,  its  addition  can  overcome  the  loss 
factor  and  reduced  degrees  of  fi-eedom.  Replication  size  must  also  be  considered  when 
taking  the  loss  factor  into  account.  For  example,  in  this  study,  internal  control  variates 
show  higher  levels  of  improvement  in  confidence  interval  width  reduction  over  the  other 
control  variates  when  rephcations  are  increased  from  10  to  20  (see  section  4.5.1). 
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However,  even  at  10  replications,  internal  control  variates  performed  at  similar  levels  of 
variance  reduction  as  the  other  control  variates.  It  is  this  problem  that  leads  to  the  process 
of  internal  control  variate  selection  discussed  above  in  section  5.2.  Analytical  and  external 
control  variates  have  fixed  loss  factors  and  associated  degrees  of  freedom  of  (n-2)/(n-3) 
and  n-2  respectively.  Despite  this  possible  disadvantage  for  internal  control  variates,  when 
some  effort  is  made  to  find  the  best  subsets,  confidence  interval  width  reduction  is 
comparable  with  that  of  the  other  control  variates  in  this  study. 

5.6  Bias  of  Analytical  Control  Variate  Point  Estimates 

Both  this  study  and  Tomick  (1988)  reveal  a  problem  with  bias  in  the  estimation  of 
system  responses  when  using  analytical  control  variates.  This  bias  is  induced  when  the 

A. 

Optimal  control  coefiBcient  is  estimated  [13].  For  the  regression  method  of  estimating 

A 

P  ,  we  assume  that  the  response  variable  and  control  variate(s)  have  a  multivariate  normal 
distribution  [12].  If  this  is  not  the  case,  bias  could  be  the  outcome  when  making  point 
estimates  of  the  response  variable. 

To  investigate  the  normality  assumption,  a  histogram  of  the  1,000  analytical  control 
variates  generated  to  estimate  system  sojourn  time  for  model  Oi  with  service  time  setting 
A.  and  transition  probability  matrix  P2  is  provided  by  figure  5.1  below.  From  the 
histogram,  it  appears  the  analytical  control  variates  may  have  more  values  to  the  right  of 
the  mean  than  might  be  expected  for  a  normal  distribution.  This  additional  weight  in  the 
higher  value  tail  could  be  the  cause  for  the  consistently  under  estimating  of  mean  sojourn 
time  by  the  analytical  control  variates.  A  chi-square  Goodness  of  Fit  test,  with  a  null 
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hypothesis  of  normality,  provides  a  test  statistic  with  a  p-value  of  approximately  .28.  This 
value  is  insufficient  to  reject  the  null  hypothesis,  however  distribution  selection  software 
indicates  other  distributions  which  provide  a  better  fit  of  the  data.  In  terms  of  the  chi- 
square  test,  the  lognormal,  chi-square,  logistic,  and  Erlang  distributions  all  provide 
significantly  better  test  statistics.  Other  design  points  provide  similar  results.  Given  an 
indication  of  non-normality  for  the  analytical  control  variates,  one  possible  means  of 
reducing  the  induced  bias  is  to  perform  a  transformation  on  the  analytical  control  variates 
that  will  produce  values  that  are  more  nearly  normally  distributed. 


Estimated  Sojourn  Time  (Analytical  Model) 


Figure  5.1.  Histogram  of  1,000  analytical  control  variate  values  of  estimated  sojourn  time, 
for  model  Qi  at  service  time  setting  A.  and  transition  probability  matrix  P2. 

Another  means  of  possibly  reducing  the  bias,  suggested  by  Dietz  [4],  is  to  adjust  the 
means  of  the  analytical  model  parameters  so  that  the  output  of  the  analytical  model  is 
closer  to  the  output  of  the  simulation  model.  Experimental  results  (see  tables  4.23-4.26) 
show  a  consistently  longer  estimate  of  system  sojourn  time  and  a  consistently  smaller 
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CPU  utilization  rate  for  the  analytical  model  versus  the  simulation  result.  In  this  study,  the 
selection  of  the  mean  for  the  exponentially  distributed  service  times  of  the  analytical  model 
is  arbitrarily  set  to  equal  the  mean  of  the  simulation  Wiebull  service  time  distributions. 
Dietz  proposes  making  the  adjustment  to  the  exponential  distributions  based  on  the 
difference  between  the  analytical  mean  queueing  length  and  an  estimate  of  the  true  average 
queueing  length.  This  method  could  be  accomplished  in  our  case  by  using  some  pilot  runs 
of  the  simulation  model  to  make  the  appropriate  adjustments.  Although  this  method  will 
produce  analytical  solutions  closer  to  the  simulation  model,  it  doesn’t  guarantee  the 
elimination  of  bias. 

Nelson  (1990)  evaluates  several  techniques  for  computing  the  controlled  response  in 
an  effort  to  remedy  bias  problems.  For  the  replication  sizes  used  in  this  study.  Nelson 
recommends  using  a  technique  known  as  splitting  for  estimating  the  response  variable  and 
associated  confidence  interval  when  bias  is  suspected.  Simply  put,  splitting  consists  of 

A. 

computing  P  fi-om  a  preliminary,  or  pilot,  sample  rather  than  the  10(20)  replications  used 
to  estimate  the  responses.  Then  even  if  the  response  variable  and  control  variate(s)  are 
not  normal,  the  controlled  estimate  will  still  be  unbiased.  To  see  this,  let  be  the  usual 
regression  estimator  of  >9  for  an  independent  sample  of  m  replications  for  response  variable 
Y  and  control  variate  vector  C.  Then  for  n  different  (independent)  replications,  using 
equation  (3.6),  the  control  variate  estimator 

F(;ff*)  =  F(«)->ff*'(C-A/,)  (5.1) 

is  unbiased  since  P  *  and  C  are  independent  of  each  other. 
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Nelson  proposes  the  use  of  what  he  calls  an  extreme  form  of  splitting  where  n 


replications  are  split  into  n  groups.  In  effect  this  consists  of  estimating  a  different  P  for 
every  y-th  replication  from  the  n-\  other  replications.  He  goes  on  to  provide 
computational  formulas  to  find  the  controlled,  split  point  estimates,  and  estimated  variance 
and  confidence  interval.  Nelson  reports  excellent  results  using  this  splitting  technique  with 
variance  estimates  nearly  identical  to  those  using  the  regression  method. 

5. 7  Recommendations  for  Future  Research 

There  are  many  possible  benefits  from  using  analytical  control  variates.  Before  we 
can  take  advantage  of  them  in  many  real  applications,  a  solution  to  the  bias  problem  must 
be  found.  Further  research  into  the  actual  distribution  of  the  analytical  control  variates 
and  associated  transformation  techniques  could  provide  significant  reduction  in  the 
observed  bias.  The  results  from  Nelson  (1990)  indicate  future  research  into  the 
application  of  the  splitting  technique  with  analytical  control  variates  should  be  pursued. 
Additionally,  the  method  of  “correcting”  the  anal3^ical  response  by  adjusting  the 
parameters  of  the  analytical  model  should  be  explored.  A  study  that  explores  the  results 
from  each  of  these  methods  separately  and  together  should  prove  to  be  very  valuable. 

If  any  of  these  techniques  eliminates  the  bias  problem,  research  into  real  world 
applications  would  be  appropriate.  One  such  example  that  could  use  the  same  mean  value 
analysis  program  ForQue  would  be  an  Air  Force  manpower  requirement  study.  Studies  of 
this  type  are  conducted  using  the  Logistics  Composites  Model  (LCOM)  computer 
simulation  model.  Jenkins  (1994)  developed  a  mean  value  heuristic  for  the  analysis  of  the 
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aircraft  sortie  generation  process,  the  same  process  used  to  evaluate  manpower 
requirements.  These  studies  showed  excellent  correlation  between  the  output  of  LCOM 
and  the  heuristic  (contained  in  ForQue)  for  a  simple  sortie  generation  process.  This 
heuristic  should  serve  as  an  excellent  analytical  model  for  generating  analytical  control 
variates  for  a  manpower  study. 

Another  area  for  future  research  is  to  find  methods  to  use  the  variance  reduction 
obtained  from  analytical  control  variates  to  validate  the  analytical  model  as  an  adequate 
surrogate  model  for  the  simulation  model.  A  technique  might  be  developed  similar  to  that 
used  to  create  response  surface  methodology  meta-models. 
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Appendix  A  :  Computer  Source  Code 


A.  1  SLAM  II  Source  Code  for  Model  Mj 

;  XX (I)  Variable  definitions 

•  'k'k'k-k’iric^’k-^-ifk-k'k-ir'k’^’ieir-iv-ir-k-k'iv-k-k-k-k-k-kifk-ir-k-k'k'k-k-^-k'k-k'k'k-ic-ic'k’k-ir'k-k'k 

r 

}  XX (1) =service  time  one  distribution  beta  parameter 
;  XX (2 ) =service  time  one  distribution  alpha  parameter 
;  XX (3) =service  time  two  distribution  beta  parameter 
;  XX (4) =service  time  two  distribution  alpha  parameter 
;  XX (5) =service  time  three  distribution  beta  parameter 
;  XX { 6) =service  time  three  distribution  alpha  parameter 
;  XX{7)=mean  of  service  time  one  distribution 

;  XX(8)=mean  of  service  time  two  distribution 

;  XX(9)=mean  of  service  time  three  distribution 

;  XX (10) =standardi2ed  work  variable  one 
;  XX (11) =standardized  work  variable  two 
;  XX ( 12 ) =standardized  work  variable  three 
;  XX (13) =standardized  arrival  variable 
;  XX(14)=sum  of  standardized  work  variable  one 

;  XX(15)=sum  of  standardized  work  variable  two 

;  XX(16)=sum  of  standardized  work  variable  three 

;  XX(17)=sum  of  standardized  arrival  variable 

;  XX (18) =arrival  time  of  previous  customer 
;  XX (25) =counter  of  customers 

;  XX (26) =standard  deviation  of  service  time  one  distribution 

;  XX (27) =standard  deviation  of  service  time  two  distribution 

;  XX (28) =standard  deviation  of  service  time  three  distribution 

f 

;  Atrib  definitions 

;  Atrib ( 1 ) =arri val  time 
;  Atrib (2) =center  one  service  time 
;  Atrib (3) =center  two  service  time 
;  Atrib (4)=center  three  service  time 
;  Atrib (5) =interarrival  time 
;  Atrib (6) =sojourn  time 

r 

GEN, I RISK, MICKEY  D  TEST  MODEL, 10/14/1995, 1000, Y,N,Y/Y,N,N/1, 132; 
LIMITS, 3, 6, 100; 

INTLC,XX(1)=6000,XX(2)=3.085,XX(3)=6000,XX(4)=3.0489; 

INTLC,XX(5)=3000,XX(6)=2.8355,XX(7)=15,XX(8)=15.5,XX(9)=15; 

INTLC,XX(10)=0,XX(11)=0,XX(12)=0,XX(13)=0,XX(14)=0,XX(15)=0; 

INTLC,XX(16)=0,XX(17)=0,XX(18)=0,XX(25)=0; 

INTLC, XX (26) =5. 3166, XX (27) =5. 5524, XX (28) =5. 7344; 

NETWORK; 

CREATE,EXPON(20, 1)  ,  ,  1,1000,  1; 

ACTIVITY; 

/ 

;  *Assign  service  times* 

ASNl  ASSIGN, ATRIB(2)=WEIBL(XX(1) ,XX(2) ,2) , ATRIB (3) =WEIBL (XX (3) , 
XX(4) ,3) ,ATRIB(4)=WEIBL (XX(5) ,XX(6) ,4) ,1; 

ACTIVITY, , ,ASN2; 
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f 

ASN2  ASSIGN, ATRIB(5)=TNOW-XX(18) ,XX(18)=TN0W, 1; 
ACTIVITY, , ,ORDR; 


ORDR 

QUEUE(l),,,; 

ACTIVITY ( 1 ) / 1 , ATRIB ( 2 ) , , PAY; 

Service 

center 

one 

r 

PAY 

QUEUE (2),,,; 

ACTIVITY(l) /2,ATRIB(3) , ,PU; 

Service 

center 

two 

} 

PU 

QUEUE (3) ,,2,BLOCK; 

ACTIVITY (1) /3,ATRIB( 4) , ,ASN3; 

Service 

center 

three 

;  *Calc  sojourn  time  and  count  customers* 

/ 

ASN3  ASSIGN, ATRIB(6)=TN0W-ATRIB{1) ,XX(25)=XX(25)  +1,1; 

ACTIVITY, , XX (25) .LE. 300, TERM;  Warm-up  period 

ACTIVITY, , ,ASN4; 

/ 

;  *Calculate  internal  cv' s* 

/ 

ASN4  ASSIGN, XX(10)=ATRIB{2)/XX(26)-XX(7)/XX(26) ,XX (11) =ATRIB (3) /XX(27) 
XX (8) /XX (27) ,XX(12)=ATRIB(4)/XX(28) -XX ( 9 ) /XX (28 ) , 
XX(13)=ATRIB(5)/20.0-1.0; 

/ 

ACTIVITY, , ,ASN5; 

} 

;  *Sum  internal  cv' s* 

r 

ASN5  ASSIGN, XX(14)=XX(14)+XX(10) ,XX(15)=XX(15)+XX(11) , 
XX(16)=XX(16)+XX(12) ,XX(17)=XX(17)+XX(13) ; 

ACTIVITY, , ,C0L1; 

/ 

;  *Collect  statistics* 

} 

COLl  COLCT(l) ,ATRIB(2) , ORDERTIME, , 1; 

ACTIVITY, , , COL2 ; 

COL2  COLCT ( 2 ) , ATRIB ( 3 ) ,  PAYTIME,  ,  1 ; 

ACTIVITY, , ,COL3; 

COL3  COLCT (3) , ATRIB (4) ,PU  TIME,,1; 

ACTIVITY, , ,COL4; 

COL4  COLCT ( 4 ) , ATRIB ( 5 ) , INTER  ARRIVE,  ,  1 ; 

ACTIVITY, , ,COL5; 

COLS  COLCT (5) ,ATRIB(6) , SOJOURN, , 1; 

ACTIVITY, , ,COL6; 

COL6  COLCT(6),XX(10),WORK  VAR1,,1; 

ACTIVITY, , ,COL7; 

COL7  COLCT(7),XX(ll),WORKVAR2,,l; 

ACTIVITY, , ,COL8; 

COL8  COLCT (8) , XX (12) , WORK  VAR3,,1; 

ACTIVITY, , ,COL9; 

COL9  COLCT ( 9 ) , XX ( 13 ) , STD  ARRV, , 1 ; 

ACTIVITY, , ,TERM; 

TERM  TERMINATE; 

END; 

FIN; 
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A.  2  FORTRAN  Subroutine  for  Model  Mi 


*  'k 

*  Main  program  * 

*  * 


PROGRAM  MAIN 
DIMENSION  NSET(5000) 

INCLUDE  'PARAM.INC 

COMMON/SCOMl/ATRIB(MATRB) ,  DD(MEQT),  DDL(MEQT),  DTNOW,  II,  MFA, 
1MST0P,NCLNR,  NCRDR,  NPRNT,  NNRUN,  NNSET,  NTAPE,  SS (MEQT) , 
2SSL(MEQT) ,TNEXT,  TNOW,  XX (MMXXV) 

COMMON  QSET{5000) 

EQUIVALENCE  (NSET ( 1) , QSET (1) ) 

NNSET=5000 

NCRDR=5 

NPRNT=6 

NTAPE=7 

OPEN (UNIT=NCRDR, FILE= ' fort . 5 ’ ) 

OPEN ( UNI T=NPRNT , FI LE= ' f ort . 6 ’ ) 

CALL  SLAM 
STOP 

END 

kkkkkkkkkk'kkk’kkkkkkkkkkkkkkkkkkkkirkkkkkkkkkkkkk'kkkkkkkkkifkkkk’k'kkk 

*  Subroutine  to  calculate  stats  and  controls  at  end  of  each  run  * 

kkkkkkkkkk'kkk'kkkirkkkkkkkkkkk'kkkkkkkkkkkk'kkkkkkkkkkkk'kkkkkkkkkkkkk 

*  Subroutine  variable  definitions 

*  Y(J)=average  time  in  system  for  J'th  replication 

*  W ( J, K) =standard  work  variable  K  for  J'th  replication 

*  S ( J, K) =service  time  mean  for  center  K  for  J'th  replication 

*  A( J) =analytical  control  variate  for  J'th  replication 

*  DENOM(J)=j 'th  denominator  for  calculation  of  external  control  variate 

*  using  Jackson  network  approximation 

*  R(J)=j'th  mean  interarrival  time 

*  Nuinber=number  of  replications 

SUBROUTINE  OTPUT 
INCLUDE  'PARAM.INC 

C0MM0N/SC0M1/ATRIB(MATRB) ,  DD(MEQT),  DDL (MEQT),  DTNOW,  II,  MFA, 
1MST0P,NCLNR,  NCRDR,  NPRNT,  NNRUN,  NNSET,  NTAPE,  SS(MEQT), 
2SSL(MEQT) ,TNEXT,  TNOW,  XX (MMXXV) 

DOUBLE  PRECISION  Y(IOOO),  W(1000,3),  S(1000,3) 

DOUBLE  PRECISION  A(IOOO),  R(IOOO) 

DOUBLE  PRECISION  DENOM(3),  SA(IOOO) 

INTEGER  matOpen,  mxCreateFull,  matClose 
INTEGER  mxGetPr,  MatPut  Matrix 
INTEGER  b,  c,  FP,  STAT 
INTEGER  J,  K,  NUMBER 
NUMBER  =  1000 
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*  Calculate  average  time  in  system  for  each  replication 

Y(NNRUN)  =  CCAVG(5) 

*  Calculate  standardized  work  variables  and  service  time  means 

DO  10  J  =  1,  3 

W(NNRUN,J)  =  XX(J+13)/SQRT(CCNUM(J+5) ) 

S (NNRUN, J)  =  CCAVG ( J) 

10  CONTINUE 

*  Calculate  mean  interarrival  time  and  standardized  arrival  variable 

SA(NNRUN)  =  XX(17)/SQRT(CCNUM(9) ) 

R (NNRUN)  =  CCAVG (4) 


*  Calculate  analytical  control  variate 

DO  20  K  =  1,  3 

DENOM(K)  =  (1.0/S (NNRUN, K) )  -  (1 . 0/R (NNRUN) ) 

20  CONTINUE 

A(NNRUN)  =  (1. 0/DENOM(1) )  +  ( 1 . O/DENOM ( 2 ) )  +  ( 1 . O/DENOM ( 3 ) ) 

*  Output  data  in  MATLAB  format 


IF  (NNRUN  ,GE.  NUMBER)  THEN 

FP  =  matOpen ( 'onethouw.mat' ,  'w') 
h-  mxCreateFull (NUMBER, 1, 0) 

call  mxCopyRealSToPtr (Y,  mxGetPr(b),  NUMBER) 
call  mxSetName(b,  'sojourn') 
stat  =  matPutMatrix (fp,  b) 

c=  mxCreateFull (NUMBER, 3, 0) 

call  mxCopyRealSToPtr (W,  mxGetPr(c),  3*NUMBER) 
call  mxSetName(c,  'stdwork') 
stat  =  matPutMatrix (fp,  c) 
call  mxFreeMatrix (c) 

call  mxCopyRealSToPtr (SA,  mxGetPr(b),  NUMBER) 
call  mxSetName(b,  'stdarrive') 
stat  =  matPutMatrix (fp,  b) 

call  mxCopyRealSToPtr  (A,  mxGetPr(b),  NUMBER) 

call  mxSetName(b,  'analytical') 

stat  =  matPutMatrix (fp,  b) 

stat  =  matClose(fp) 

call  mxFreeMatrix (b) 

END  IF 

RETURN 

END 
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A.  3  SLAM II  Source  Code  for  Model  M3 


;  XX (I)  Variable  definitions 

•  ^•k-k^-k-k-k'k-k'k-k-k-k-k^-k^^-k'k-k-ir-if-k-k-k-h-k-k-k-k-k-k-k-k-k-k-k-k^-k-k-k-k-k-k-k-k-kit'k 

;  XX ( 1) =exponential  service  time  one  distribution  mean 
;  XX (2)=exponential  service  time  two  distribution  mean 
;  XX (3) =exponential  service  time  three  distribution  mean 

•  'k'k'k'k'k‘k‘^'k'k-k'k'k’k’k-k’k'k‘k'if'^f'k'k'k’k'k’ir'-)e-if)ifk'k’k’k'ie-k-k'k-k-k-k‘k'k'k'k‘k'k‘k‘k-k'ie'k 

;  Atrib  definitions 

;  Atrib (1) =arrival  time 
;  Atrib (2 ) =center  one  service  time 
;  Atrib (3) =center  two  service  time 
;  Atrib (4) =center  three  service  time 
;  Atrib (5) =sojourn  time 

GEN, I RISK, MICKEY  D  TEST  MODEL,10/14/1995,1000,Y,N,y/Y,N,N/l,132; 
LIMITS, 3, 5, 200; 

INTLC,XX(1)=15.0,XX{2)=15.5,XX(3)=15.0; 

INTLC,XX{25)=0; 

NETWORK; 

/ 

CREATE, EXPON (20,1)  ,  ,1,  1000,1; 

ACTIVITY; 

/ 

;  *Assign  service  times* 

/ 

ASNl  ASSIGN,ATRIB(2)=EXP0N{XX(1) ,2) ,ATRIB{3)=EXPON(XX(2) ,3) , 
ATRIB(4)=EXPON(XX{3) ,4) ; 


ACTIVITY, , , ORDR; 

/ 

ORDR 

QUEUE (1),,,; 

ACTIVITY (1) /I, ATRIB (2) ,,PAY; 

Service 

center 

one 

/ 

PAY 

QUEUE(2),,,; 

ACTIVITY (l)/2, ATRIB (3) ,,PU; 

Service 

center 

two 

PU 

QUEUE (3),,,; 

ACTIVITY(l) /3,ATRIB(4) ,  ,ASN2; 

Service 

center 

three 

;  *  Calculate  sojourn  time  and  count  customers* 

} 

ASN2  ASSIGN,ATRIB(5)=TNOW-ATRIB(l) ,XX(25)=XX(25)  +1; 

ACTIVITY, , XX (25) .LE. 3 00, TERM;  Warm-up  period 

ACTIVITY, , XX (25) . GT . 300, COLl; 

f 

COLl  COLCT (1) ,  ATRIB (5) , SOJOURN, , 1; 

ACTIVITY, , , TERM; 

/ 

TERM  TERMINATE; 

END; 

FIN; 
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A.  4  FORTRAN  Subroutine  for  Model  M3 


*  Main  program  * 


PROGRAM  MAIN 
DIMENSION  NSET(5000) 

INCLUDE  'PARAM.INC' 

C0MM0N/SC0M1/ATRIB{MATRB) ,  DD(MEQT),  DDL(MEQT),  DTNOW,  II,  MFA, 
1MST0P,NCLNR,  NCRDR,  NPRNT,  NNRUN,  NNSET,  NTAPE,  SS{MEQT), 
2SSL(MEQT) ,TNEXT,  TNOW,  XX (MMXXV) 

COMMON  QSET(5000) 

EQUIVALENCE  (NSET (1) ,QSET (1) ) 

NNSET=5000 

NCRDR=5 

NPRNT=6 

NTAPE=7 

OPEN (UNIT=NCRDR, FILE= ' fort . 5 ’ ) 

OPEN(UNIT=NPRNT,FILE='fort.6' ) 

CALL  SLAM 

STOP 

END 

*  Subroutine  to  calculate  stats  and  controls  at  end  of  each  run  * 


*  Subroutine  variable  definitions 

*  Y(J)=average  time  in  system  for  J'th  replication 

*  Nuinber=number  of  replications 

SUBROUTINE  OTPUT 
INCLUDE  'PARAM.INC 

COMMON/SCOMl/ATRIB(MATRB) ,  DD(MEQT),  DDL(MEQT],  DTNOW,  II,  MFA, 
1MST0P,NCLNR,  NCRDR,  NPRNT,  NNRUN,  NNSET,  NTAPE,  SS(MEQT), 
2SSL(MEQT) ,TNEXT,  TNOW,  XX (MMXXV) 

DOUBLE  PRECISION  Y(IOOO) 

INTEGER  matOpen,  itixCreateFull,  matClose 
INTEGER  mxGetPr,  MatPut  Matrix 
INTEGER  b,  FP,  STAT 
INTEGER  NUMBER 

NUMBER  =  1000 


*  Calculate  average  time  in  system  for  each  replication 
Y (NNRUN)  =  CCAVG(l) 
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*  Output  data  in  MATLAB  format 

IF  (NNRUN  .GE.  NUMBER)  THEN 

FP  =  matOpen ( ' extthou.mat ' ,  'w') 

b=  mxCreateFull (NUMBER, 1, 0) 

call  mxCopyReal8ToPtr (Y,  mxGetPr(b),  NUMBER) 

call  mxSetName (b,  ' extcontrol ' ) 

stat  =  matPutMatrix (fp,  b) 

stat  =  matClose(fp) 

call  mxFreeMatrix (b) 


END  IF 
RETURN 
END 
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A.S  SLAM  II  Source  Code  for  Models  Qi  and  Qc 


XX (I)  Variable  definitions 
******■*•****  +  +  *■*****■*■■*•**********+*•*■  ++*******  +  ■****  +  +  * 

XX{l)=event  counter 

XX(2)=time  when  first  statistics  are  gathered 

XX (3) =standardized  work  variable  one 

XX(4)=suin  of  standardized  work  variable  one 

XX (5) =standardized  work  variable  two 

XX(6)=suin  of  standardized  work  variable  two 

XX (7) =standardized  work  variable  three 

XX(8)=siim  of  standardized  work  variable  three 

XX (9) =standardized  work  variable  four 

XX(10)=sum  of  standardized  work  variable  four 

XX (11) =standardized  work  variable  five 

XX(12)=sum  of  standardized  work  variable  five 

XX (13) =standardized  work  variable  six 

XX(14)=sum  of  standardized  work  variable  six 

XX (15) =customer  capacity  for  subnetwork  (=25  for  Q-1) 

XX (20)=  service  time  one  distribution  beta  parameter 

XX (21)=  service  time  one  distribution  alpha  parameter 

XX (22)=  service  time  one  distribution  mean 

XX (23)=  service  time  one  distribution  standard  deviation 

XX (24)=  service  time  three  distribution  beta  parameter 

XX (25)=  service  time  three  distribution  alpha  parameter 

XX (26)=  service  time  three  distribution  mean 

XX (27)=  service  time  three  distribution  standard  deviation 

XX (28)=  service  time  four  distribution  beta  parameter 

XX (29)=  service  time  four  distribution  alpha  parameter 

XX (30)=  service  time  four  distribution  mean 

XX (31)=  service  time  four  distribution  standard  deviation 

XX (32)=  service  time  five  distribution  beta  parameter 

XX (33)=  service  time  five  distribution  alpha  parameter 

XX (34)=  service  time  five  distribution  mean 

XX (35)=  service  time  five  distribution  standard  deviation 

XX (36)=  service  time  six  distribution  beta  parameter 

XX (37)=  service  time  six  distribution  alpha  parameter 

XX (38)=  service  time  six  distribution  mean 

XX (39)=  service  time  six  distribution  standard  deviation 

Atrib  definitions 

Atrib (1) =time  leaving  center  one 
Atrib (2) =center  two  service  time 
Atrib (3) =center  three  service  time 
ATRIB(4)=center  four  service  time 
Atrib (5) =center  five  service  time 
Atrib (6)=center  six  service  time 
Atrib (7) =center  one  service  time 
Atrib (8)=sojourn  time 

Atrib (9) =routing  random  draw  from  center  two 
Atrib (10) =initial  starting  position  random  draw 


GEN, THOMAS  H.  IRISH, COMPUTER,  11/28/1995, 1000, Y, N, Y/y,N, N/1, 132 
LIMITS, 8, 10,25; 


INTLC,XX(1)=0,XX(2)=0,XX{3)=0,XX(4)=0,XX(5)=0,XX(6)=0,XX(7)=0,XX(8)=0; 

INTLC,XX{9)=0,XX(10)=0,XX(11)=0,XX{12)=0,XX(13)=0,XX(14)=0,XX(15)=25; 

INTLC, XX (20) =1000.0, XX (21) =1.4 6824, XX (22) =100.0, XX (23) =69. 2516; 

INTLC,XX(24)=10.0,XX(25)=5.6476,XX(26)=1.39,XX(27)=.28477; 

INTLC,XX(28)=10.0,XX(29)=5.6476,XX(30)=1.39,XX(31)=.28477; 

INTLC,XX(32)=1000.0,XX(33)=2.61249,XX(34)=12.5,XX(35)=5.14217; 

INTLC,XX(36)=1000.0,XX(37)=2.61249,XX(38)=12.5,XX(39)=5.14217; 

NETWORK; 

RESOURCE/1,  SPACE (XX (15) ) ,7,8; 

r 

CREATE, .01,,,25,1; 

ACTIVITY; 

} 

;  *Send  customers  to  initial  service  centers* 


STRT  ASSIGN,ATRIB(10)=UNFRM(0,1,8),1; 

ACTIVITY, ,ATRIB( 10) .LE. .708,ASNA; 

ACTIVITY, ,ATRIB(10) .GT. . 708 -AND.NNRSC (1) .LT.1,CAP1; 

ACTIVITY, ,ATRIB(10) . GT . .708.AND-NNRSC (1) .GT. 0,CAP2; 

/ 

CAP2  AWAIT (7 ), SPACE,, 1;  (Center  2  in  Q-c) 

ACTIVITY,  ,'ATRIB(10)  .GT.  .  708  .  AND.  ATRIB  ( 10 )  .LE.  .  877, TWOS; 
ACTIVITY,, ATRIB( 10) . GT . . 877 .AND. ATRIB (10) .LE. .9078,SRV4; 
ACTIVITY, ,ATRIB(10) . GT .. 9078 .AND. ATRIB ( 10) .LE. . 9386,  SRV4; 
ACTIVITY, ,ATRIB(10) .GT. . 9386 .AND, ATRIB ( 10) .LE. .9693,SRV5; 
ACTIVITY, ,ATRIB(10) .GT. .9693,SRV6; 

9 

ASNA  ASSIGN,ATRIB(7)=WEIBL(XX(20),XX(21),1),1;  Service  time  1 
ACTIVITY, , ,Q1; 

Q1  QUEUE(l),,,; 

ACTIVITY (25) /I, ATRIB (7)  ; 

ASNl  ASSIGN, XX ( 1 ) =XX ( 1 ) +1 , ATRIB ( 1 ) =TNOW, 1 ; 

ACTIVITY, ,XX(1) ,GE. 2000, TERM; 

ACTIVITY, ,XX(1) .EQ.500; 

ACTIVITY, , XX (1)  .GT.500,WRK1; 

ACTIVITY, , XX (1) .LT.500,CAP1; 

CLKl  ASSIGN,XX(2)=TNOW, 1; 

ACTIVITY; 

9 

;  *Calc  std  work  cv  and  service  time  #1* 

9 

WRKl  ASSIGN,XX(3)=ATRIB(7)/XX(23)-XX(22)/XX(23) ,XX(4)=XX(4)  +  XX (3), 1 
ACTIVITY; 

WKCl  COLCT(7) ,XX(4) ,WORK  ONE,,l; 

ACTIVITY; 

SVCl  COLCT ( 13 ),  ATRIB (7) , SERVICE  1,,1; 

ACTIVITY, , ,CAP1; 

9 

TWOS  ASSIGN,ATRIB(2)=EXPON(l,2) ,1;  Service  time  2(3) 

ACTIVITY; 

Q2  QUEUE(2) , , , ;  Service  center  2(3) 

ACTIVITY (1) /2,ATRIB (2) ; 


Service  center  1 

Count  customers 
End  replication 
Start  stats 
Collect  stats 
Warm-up  period 
Record  stat  start  time 


;  *Draw  routing  rv,  and  count  event* 


ASN2  ASSIGN,XX(1)=XX(1)+1,ATRIB(9)=UNFRM(0,1,10),1; 

ACTIVITY, ,XX(1) .GE. 2000, TERM;  End  replication 
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ACTIVITY, ,XX(1) .EQ.500; 

ACTIVITY, , XX (1) .GT.500,WRK2; 

ACTIVITY, , XX (1) .LT.500,PROB; 

CLK2  ASSIGN,XX(2)=TN0W, 1; 

ACTIVITY; 

r 

;  *Calc  std  work  cv  and  service  time  #2(3)* 

} 

WRK2  ASSIGN,XX(5)=ATRIB(2)-1,XX(6)=XX(6)  + 

ACTIVITY; 

WKC2  C0LCT(8) ,XX(6) , WORK  TWO,, 1; 

ACTIVITY; 

SVC2  COLCT { 14 ) ,ATRIB (2) , SERVICE  2,,1; 

ACTIVITY, , , TWOS ; 

r 

;  *Branch  to  next  service  center* 

t 

PROB  G00N,1; 

ACTIVITY, ,ATRIB( 9) .LE. .2  .AND.  XX(1).LT.500; 

ACTIVITY, ,ATRIB (9) .LE. .2  .AND.  XX ( 1 ) . GE . 500 , RTEl ; 

ACTIVITY, ,ATRIB (9) .GT. .2  .AND.  ATRIB(9).LE..56,SRV3; 

ACTIVITY, ,ATRIB (9) .GT. .56  .AND.  ATRIB ( 9 ) . LE. . 92, SRV4 ; 

ACTIVITY,, ATRIB( 9) .GT. .92  .AND.  ATRIB ( 9) . LE . . 96, SRV5; 

ACTIVITY , , ATRIB ( 9 )  . GT . . 9 6 ,  SRV6 ; 

FREl  FREE , S  PACE , 1 ; 

ACTIVITY, , ,ASNA; 

RTEl  COLCT ( 19) , ALL, ROUTE  1,,1;  Collect  routing  #1  info 

ACTIVITY; 

FRE2  FREE , S  PACE , 1 ; 

ACTIVITY, , ,ASN7; 

/ 

;  *Assign  service  times  for  centers  4-6  (5-7)  and  calc  routing  info* 

/ 

SRV3  ASSIGN, ATRIB(3)=WEIBL(XX(24) ,XX(25) ,3) ,1; 

ACTIVITY, , XX (1) .LT.500,Q3; 

ACTIVITY; 

RTE3  COLCT (20) , ALL, ROUTES,  ,  1; 

ACTIVITY, , ,Q3; 

SRV4  ASSIGN,ATRIB(4)=WEIBL{XX(28) ,XX(29) ,4),1; 

ACTIVITY, , XX ( 1 ) .LT.500,Q4; 

ACTIVITY; 

RTE4  COLCT (21) , ALL, ROUTE  4,, 1; 

ACTIVITY, , , Q4 ; 

SRV5  ASSIGN,ATRIB(5)=WEIBL(XX(32) ,XX(33) ,5),1; 

ACTIVITY, , XX (1) .LT.500,Q5; 

ACTIVITY; 

RTE5  COLCT (22 ), ALL, ROUTE  5,,  1; 

ACTIVITY, , ,Q5; 

SRV6  ASSIGN,ATRIB(6)=WEIBL(XX(36) ,XX(37) ,6) ,1; 

ACTIVITY, ,XX(1) .LT.500,Q6; 

ACTIVITY; 

RTE6  COLCT (23) , ALL, ROUTE  6,,1; 

ACTIVITY, , ,Q6; 

} 

Q3  QUEUE (3),,,;  Service  center  3(4) 

ACTIVITY (1) /3, ATRIB (3) ; 

ASN3  ASSIGN,XX(1)=XX(1)+1,1;  Count  event 

ACTIVITY, ,XX(1) .GE. 2000, TERM;  End  replication 


Start  stats 
Collect  stats 
Warm-up  period 
Record  stat  start  time 


XX (5) , 1; 
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ACTIVITY, ,XX(1) .EQ.500; 
ACTIVITY, ,XX{1)  .GT.500,WRK3; 

ACTIVITY, , XX (1) .LT. 500, TWOS; 
CLK3  ASSIGN,XX(2)=TNOW, 1; 

ACTIVITY; 


Start  stats 
Collect  stats 
Warm-up  period 
Record  stat  start  time 


;  *Collect  std  work  and  service  time  #3(4) 


/ 

WRK3  ASSIGN, XX(7)=ATRIB(3)/XX{27)-XX{26)/XX(27),XX(8)=XX(8)  +  XX(7),1; 
ACTIVITY; 

WKC3  COLCT ( 9 ) , XX ( 8 ) , WORK  THREE, , 1 ; 

ACTIVITY; 

SVC3  COLCT ( 15) ,ATRIB (3) , SERVICE  3,,1; 

ACTIVITY, , , TWOS ; 

r 

Q4  QUEUE (4),,,; 

ACTIVITY ( 1) /4,ATRIB (4)  ; 

ASN4  ASSIGN,XX(1)=XX{1)+1; 

ACTIVITY, , XX (1) .GE. 2000, TERM; 

ACTIVITY,,XX(1) .EQ.500; 

ACTIVITY, , XX (1)  .GT.500,WRK4; 

ACTIVITY, , XX ( 1 ) .LT. 500, TWOS; 

CLK4  ASSIGN,XX(2)=TNOW, 1; 

ACTIVITY; 

/ 

;  *Calc  and  collect  std  work  and  service  time  #4(5)* 

WRK4  ASSIGN,XX(9)=ATRIB(4)/XX(31)-XX(30)/XX(31) ,XX(10)=XX(10)+XX(9) ,1; 
ACTIVITY; 

WKC4  COLCT (10) , XX (10) , WORK  FOUR, ,1; 

ACTIVITY; 

SVC4  COLCT ( 16) ,ATRIB (4) , SERVICE  4,,1; 

ACTIVITY, , ,TWOS; 

r 

Q5  QUEUE (5),,,; 

ACTIVITY(l) /5,ATRIB(5)  ; 

ASN5  ASSIGN,XX(1)=XX{1)+1; 

ACTIVITY, ,XX(1) .GE. 2000, TERM; 

ACTIVITY, ,XX{1) .EQ.500; 

ACTIVITY, , XX (1)  .GT.500,WRK5; 

ACTIVITY, ,XX{1) .LT. 500, TWOS; 

CLK5  ASSIGN, XX(2)=TNOW,l; 

ACTIVITY; 

} 

;  *Calc  and  collect  std  work  and  service  time  #5(6)* 

r 

WRK5  ASSIGN,XX(11)=ATRIB(5)/XX(35)-XX(34)/XX(35) , 

XX(12)=XX(12)+XX(11)  ,1; 

ACTIVITY; 

WKC5  COLCT(ll) ,XX(12) ,WORK  FIVE, ,1; 

ACTIVITY; 

SVC5  COLCT (17) ,ATRIB( 5) , SERVICE  5,,1; 

ACTIVITY, , ,TWOS; 

} 

Q6  QUEUE (6),,,;  Service  center  6(7) 

ACTIVITY(l) /6,ATRIB(6)  ; 

ASN6  ASSIGN,XX(1)=XX(1)+1;  Count  event 

ACTIVITY, , XX (1) .GE. 2000, TERM;  End  replication 


Service  center  5(6) 

Count  event 
End  replication 
Start  stats 
Collect  stats 
Warm-up  period 
Record  stat  start  time 


Service  center  4(5) 

Count  event 
End  replication 
Start  stats 
Collect  stats 
Warm-up  period 
Record  stats  start  time 
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ACTIVITY, , XX (1) .EQ.500; 
ACTIVITY, ,XX(1)  .GT.500,WRK6; 

ACTIVITY, , XX ( 1 ) . LT . 5  0  0 , TWOS ; 
CLK6  ASSIGN,XX(2)=TNOW, 1; 

ACTIVITY; 


Start  stats 
Collect  stats 
Warm-up  period 
Record  stat  start  time 


;  *Calc  and  collect  std  work  and  service  time  #6(7)* 


WRK6 

ASSIGN,XX(13)=ATRIB(6)/XX(39)-XX(38)/XX(39) , 
XX ( 14 ) =XX ( 14 ) +XX ( 13) ,  1 ; 

ACTIVITY; 

WKC6 

C0LCT(12) ,XX(14) ,WORK  SIX, ,1; 

ACTIVITY; 

SVC  6 

C0LCT(18) ,ATRIB(6) , SERVICE  6,,1; 

ACTIVITY, , , TWOS ; 

} 

ASN7 

AS  S I GN , ATRI B ( 8 ) =TNOW-ATRI B ( 1 ) , 1 ; 

ACTIVITY; 

COLl 

COLCT ( 1 ) , ATRIB ( 8 ) , RESPONSE  TIME, ,  1 ; 

Sojourn  time 

ACTIVITY,  ,  ,ASNA; 

r 

TERM 

TERMINATE,  1; 

} 

CAPl 

FIN; 

AWAIT (8) , SPACE, ,1;  (Service 

ACTIVITY, , ,TWOS; 

END; 

center  2  in  Q-c) 
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A.6  FORTRAN  Subroutine  for  Models  Qi  and  Qc 


*  Main  program  * 


PROGRAM  MAIN 
DIMENSION  NSET{5000) 

INCLUDE  'PARAM.INC 

COMMON/SCOMl/ATRIB(MATRB) ,  DD(MEQT),  DDL(MEQT),  DTNOW,  II,  MFA, 
1MST0P,NCLNR,  NCRDR,  NPRNT,  NNRUN,  NNSET,  NTAPE,  SS(MEQT), 
2SSL{MEQT) ,TNEXT,  TNOW,  XX(MMXXV) 

COMMON  QSET{5000) 

EQUIVALENCE  (NSET (1)  ,  QSET (1) ) 

NNSET=5000 

NCRDR=5 

NPRNT=6 

NTAPE=7 

OPEN (UNIT=NCRDR, FILE= ' fort . 5 ' ) 

OPEN(UNIT=NPRNT,FILE='fort.6’  ) 

OPEN (10,  FILE= 'analytical. in ’ ,STATUS=' UNKNOWN ' ) 

CALL  SLAM 
CLOSE (10) 

STOP 

END 

*  Subroutine  to  calculate  stats  and  controls  at  end  of  each  run  * 
*■***************•*■***+**********•*******•************************+** 

*  Subroutine  variable  definitions 

*  Y(NNRUN)  =  average  return  time  for  each  replication  (response) 

*  WAIT (NNRUN, J)  =  average  wait  at  J'th  station  for  each  replication 

*  UTILIZ (NNRUN)  =  average  utilization  of  station  2  (CPU) 

*  W (NNRUN, J)  =  Standardized  work  variable  for  J'th  station 

*  R (NNRUN, J)  =  standardized  routing  variable  for  J'th  route 

*  PROB(J, K)  =  empirical  transition  probability  matrix 

*  S (NNRUN, J)  =  average  service  time  for  J'th  station 

*  DENOM(J)  =  denominator  for  J'th  standardized  routing  variable 

*  P(J)  =  actual  routing  probability 

*  NUMBER  =  number  of  replications 

*  CUST  =  number  of  customers  in  network 

*  STATION  =  number  of  stations  in  network 

SUBROUTINE  OTPUT 
INCLUDE  'PARAM.INC 

COMMON/SCOMl/ATRIB(MATRB) ,  DD(MEQT),  DDL(MEQT),  DTNOW,  II,  MFA, 
1MST0P,NCLNR,  NCRDR,  NPRNT,  NNRUN,  NNSET,  NTAPE,  SS(MEQT), 
2SSL(MEQT) ,TNEXT,  TNOW,  XX(MMXXV) 

DOUBLE  PRECISION  ¥(1000),  W(1000,6),  R(1000,5),  WAIT(1000,5) 
DOUBLE  PRECISION  UTILIZ (1000) ,  DENOM(5),  P(5) 

REAL  PROB(6,6),  S(1000,6) 

INTEGER  matOpen,  mxCreateFull,  matClose 
INTEGER  mxGetPr,  MatPut  Matrix 
INTEGER  b,  c,  d,  FP,  STAT 
INTEGER  J,  K,  NUMBER,  CH(6) 

INTEGER  CUST,  STATION 
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DATA 

CH  /  25,1,1, 

1,1,1 

/ 

DATA 

(PROBd,  J)  , 

J=l,6) 

/ 

0.0,  1.0,  0.0,  0.0, 

0. 

0, 

+ 

(PROB(2, J) , 

J=l,6) 

/ 

0.,  0.,  0.,  0.,  0., 

0. 

/ 

+ 

(PROB(3, J) , 

J=l,6) 

/ 

0.,  1.,  0.,  0.,  0., 

0. 

/ 

+ 

(PROB(4, J) , 

J=l,6) 

/ 

0.,  1.,  0.,  0.,  0., 

0. 

/ 

+ 

(PROB(5, J) , 

J=l,6) 

/ 

0.,  1.,  0.,  0.,  0., 

0. 

/ 

+ 

(PROB(6, J) , 

J=l,6) 

/ 

0. ,  1. ,  0. ,  0. ,  0. , 

0. 

/ 

P(l)  =  .2 
P(2)  =  .36 
P(3)  =  .36 
P(4)  =  .04 
P{5)  =  .04 

NUMBER  =  1000 
GUST  =25 
STATION  =  6 

*  Calculate  system  response  characteristics  * 


Y(NNRUN)  =  CCAVG(l) 

UTILIZ (NNRUN)  =  (CCAVG ( 14 ) *CCNUM( 14 ) ) / (TNOW-XX (2 ) ) 

*  Calculate  internal  control  variates  * 


*  Calculate  standardized  work  variables 

DO  10  J  =  1,  6 

W(NNRUN,J)  =  XX{2*J+2)/SQRT(CCNUM(J+6) ) 

10  CONTINUE 

*  Calculate  standardized  routing  variables 

DENOM(l)  =  SQRT(CCNUM(8) * (l-P(l) ) *P(1) ) 

R(NNRUN,1)  =  (CCNUM(7)-CCNUM(8)*P(1) )/DENOM{l) 

DO  20  J=2,  5 

DENOM{J)  =  SQRT(CCNUM(8)*(1-P(J))*P(J) ) 

R (NNRUN, J)  =  (CCNUM(J+7)-CCNUM(8) *P (J) ) /DENOM{J) 
20  CONTINUE 

*  Output  response  data  and  internal  cv's  in  MATLAB  format 

IF  (NNRUN  .GE.  NUMBER)  THEN 

FP  =  matOpen ( 'model. mat' ,  'w') 

b=  rnxCreateFull  (NUMBER,  1,  0) 

call  mxCopyReal8ToPtr (Y,  mxGetPr(b),  NUMBER) 
call  mxSetName(b,  'sojourn') 
stat  =  matPutMatrix (fp,  b) 

call  mxCopyReal8ToPtr (UTILIZ,  mxGetPr(b),  NUMBER) 
call  mxSetName(b,  'utiliz') 
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stat  =  matPutMatrix (fp,  b) 
call  mxFreeMatrix (b) 

c=  mxCreateFull (NUMBER, 6, 0) 

call  mxCopyRealSToPtr (W,  mxGetPric),  6*NUMBER) 
call  mxSetName (c,  'stdwork') 
stat  =  matPutMatrix (fp,  c) 
call  mxFreeMatrix (c) 

d=  mxCreateFull (NUMBER, 5,0) 

call  mxCopyRealSToPtr (R,  mxGetPr(d),  5*NUMBER) 
call  mxSetName (d,  'stdroute') 
stat  =  matPutMatrix (fp,  d) 

call  mxCopyReal8ToPtr (WAIT,  mxGetPr(d),  5*NUMBER) 
call  mxSetName (d,  'wait') 
stat  =  matPutMatrix (fp,  d) 
call  mxFreeMatrix (d) 

stat  =  matClose(fp) 


END  IF 

*  Calculate  and  format  output  for  analytical  control  program  (ForQue)  * 
**************■**********•*•*•**********•*•■*•*************•*********  +  *********** 

*  Calculate  average  service  times  and  actual  routing  proportions 

DO  30  J=l,  6 

S(NNRUN,J)  =  CCAVG(J+12) 

30  CONTINUE 
K  =  2 

PROB(K,l)  =  CCNUM(19)/CCNUM(8) 

DO  40  J=3,  6 

PROB(K, J)  =  CCNUM(J+17) /CCNUM(8) 

40  CONTINUE 

*  Output  data  for  analytical  program  for  each  replication 

WRITE  (10,100)  OUST,  STATION 
WRITE  (10,200) 

WRITE  (10,500)  (S(NNRUN,J),  J=l,6) 

WRITE  (10,300)  (CH(J),  J=l,6) 

WRITE  (10,200) 

WRITE  (10,400)  ((PROB(J,K),  K=l,6),  J=l,6) 

WRITE  (10,200) 

100  FORMAT  (213) 

200  FORMAT  {/) 

300  FORMAT  (613) 

400  FORMAT  (6(6F9.5/)) 

500  FORMAT  (6F9.5) 

RETURN 

END 
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A.  7  SLAM  II  Source  Code  for  Model  Q2 


•  •k-k'k'k'k-k'kifk-k-k'k’k’k'k-k'k-k-k-k’k-k’k'k’k’k’k-if'k'ie'k-k'k’k'ick'k'k-k^’k-k-k-k-k-k-ifk-k'k'k 

;  XX (I)  Variable  definitions 

•  •k’k-k-k'k'k-k'k'k^-k-k’k’k-ir'k'icie'ifk'k'k-k-k-k-k'k'k-k-k-k-k-k^-k-k’k'k’k'k'k'k-k'lr-k-k’k’k-k-k-k 

;  XX(l)=event  counter 

;  XX{2)=time  when  first  statistics  are  gathered 
;  Atrib  definitions 

•  ’k'k-k-k-ir'if’kifk'k'k'k'it’k'k'k'k-kir-kir-ie'^-k-k-k-k-k-k-k-k-k-ir-k'k'k-k'k'fflf'kir'lr'k'k’k’ifk'k'k'k 

;  Atrib (1) =time  leaving  center  one 
;  Atrib (2) =center  two  service  time 
;  Atrib (3) =center  three  service  time 
;  Atrib ( 4 ) =center  four  service  time 
;  ATRIB (5) =center  five  service  time 
;  ATRIB ( 6) =center  six  service  time 
;  ATRIB (7) =center  one  service  time 
;  ATRIB(8)=sojourn  time 

;  ATRIB ( 9) =routing  random  draw  from  center  two 
;  ATRIB (10) =initial  starting  position  random  draw 


GEN, THOMAS  H.  IRISH,  COMPUTER,  11/28/1995, 1000, Y, N, Y/Y, N, N/ 1 , 132 ; 
LIMITS, 6, 10, 25; 

INTLC,XX(1)=0,XX{2)=0; 

NETWORK; 

} 

CREATE,0,, ,25,1; 

ACTIVITY; 


;  *  Send  customers  to  initial  service  center* 


STRT  ASSIGN,ATRIB(10)=UNFRM(0,1,8) ,1; 

ACTIVITY,  ,ATRIB(10)  .LE. ,708,ASNA; 

ACTIVITY, ,  ATRIB (10)  .GT. . 708 .AND.  ATRIB ( 10 )  .LE. .877,TWOS; 
ACTIVITY, , ATRIB (10)  .GT. . 877 .AND.  ATRIB ( 10) .LE. .9078,SRV4; 
ACTIVITY, ,ATRIB(10) .GT. . 9078 .AND. ATRIB ( 10) .LE. . 9386,SRV4; 
ACTIVITY, , ATRIB (10) .GT. . 9386. AND. ATRIB (10) .LE. .9693,SRV5; 
ACTIVITY, , ATRIB (10) .GT. .9693,SRV6; 


ASNA  ASSIGN, ATRIB (7) =EXP0N (100, 1),1;  Service  time  1 

ACTIVITY, , ,Q1; 


Q1  QUEUE (1),,,; 

ACTIVITY (25) /I, ATRIB (7) ; 

ASNl  ASSIGN, XX ( 1 ) =XX ( 1 ) +1 ,  ATRIB ( 1 ) =TNOW, 1 ; 
ACTIVITY, ,XX(1) .GE.2000,TERM; 
ACTIVITY, ,XX(1) .EQ.500; 
ACTIVITY,,XX(1)  .GT.500,ZAAH; 
ACTIVITY, , XX (1) .LT. 500, TWOS; 

CLKl  ASSIGN,XX(2)=TNOW,l; 

ACTIVITY; 

ZAAH  G00N,1; 

ACTIVITY, , ,TWOS; 

r 

TWOS  ASSIGN,ATRIB(2)=EXPON(l,2) , 1; 
ACTIVITY; 


Service  center  1 

Count  event 
End  replication 
Start  stats 
Collect  stats 
Warm-up  period 
Record  stat  start  time 


Service  time  2 
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Q2 


QUEUE (2),,,; 

ACTIVITY (1) /2,ATRIB{2) ; 


Service  center  2 


*Count  event  and  assign  branching  rv* 


ASN2  ASSIGN,XX{1)=XX(1) +1, ATRIB (9) =UNFRM(0, 1, 10)  ,  1; 


CLK2 


SVC2 


ACTIVITY, , XX ( 1 )  . GE. 2000,  TERM; 
ACTIVITY, ,XX(1) .EQ.500; 

ACTIVITY, ,XX(1) .GT.500,SVC2; 
ACTIVITY, , XX ( 1 ) .LT.500, PROB; 

AS  S I GN , XX { 2 ) =TNOW , 1 ; 

ACTIVITY; 

COLCT(14) , ATRIB (2) , SERVICE  2,,1; 
ACTIVITY, , ,TWOS; 


End  replication 
Start  stats 
Collect  stats 
Warm-up  period 
Record  stat  start 


time 


PROB  GOON,l; 

;*  Branch  customers  to  service  center  3-6* 
f 

ACTIVITY,  , ATRIB (9)  .LE. .2  .AND.  XX(1).LT.500,ASNA; 
ACTIVITY, , ATRIB (9) .LE. .2  .AND.  XX ( 1 ) . GE. 500; 

ACTIVITY, , ATRIB (9) .GT. .2  .AND.  ATRIB ( 9) . LE. . 56, SRV3 ; 
ACTIVITY, , ATRIB (9)  .GT. .56  .AND.  ATRIB(9)  .LE. . 92,  SRV4; 
ACTIVITY, , ATRIB (9)  .GT. .92  .AND.  ATRIB ( 9)  . LE. . 96,  SRV5; 
ACTIVITY, , ATRIB (9)  .GT. .96,SRV6; 


;  *Assign  service  times* 

r 

SRV3  ASSIGN, ATRIB (3) =EXPON (1.39, 3) ,1; 
ACTIVITY, ,  ,Q3; 

SRV4  ASSIGN,ATRIB(7)=EXPON (1.39,4) ,1; 
ACTIVITY, ,  ,Q4; 

SRV5  ASSIGN,ATRIB(3)=EXPON(12.5,5) , 1; 
ACTIVITY, , ,Q5; 

SRV6  ASSIGN,ATRIB(6)=EXPON(12.5,6) ,1; 

ACTIVITY, ,  ,Q6; 
f 

Q3  QUEUE (3),,,; 

ACTIVITY(l) /3,ATRIB(3) ; 

ASN3  ASSIGN,XX(1)=XX(1)+1,1; 

ACTIVITY, ,XX (1) .GE. 2000, TERM; 
ACTIVITY, , XX (1) .EQ.500; 

ACTIVITY, ,XX(1)  .GT.500,ZAAJ; 

ACTIVITY, , XX ( 1 ) . LT . 5  0  0 , TWOS ; 

CLK3  ASSIGN,XX(2)=TNOW, 1; 

ACTIVITY; 

ZAAJ  GOON,l; 

ACTIVITY, , ,TWOS; 

f 

Q4  QUEUE (4),,,; 

ACTIVITY(l) /4,ATRIB(7) ; 

ASN4  ASSIGN, XX (1)=XX(1) +1,1; 

ACTIVITY, ,XX(1) .GE. 2000, TERM; 
ACTIVITY, , XX (1) .EQ.500; 

ACTIVITY, ,XX(1)  .GT.500,ZAAK; 

ACTIVITY, , XX (1) . LT.500, TWOS; 

CLK4  ASSIGN,XX(2)=TNOW,l; 

ACTIVITY; 


Service  center  3 

Count  event 
End  replication 
Start  stats 
Collect  stats 
Warm-up  period 
Record  stat  start  time 


Service  center  4 

Count  event 
End  replication 
Start  stats 
Collect  stats 
Warm-up  period 
Record  stat  start  time 
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ZAAK  G00N,1; 

ACTIVITY, , , TWOS ; 

} 

Q5  QUEUE (5),,,; 

ACTIVITY(l) /5,ATRIB(5) ; 

ASN5  ASSIGN, XX (1)=XX(1) +1,1; 

ACTIVITY,, XX (1) .GE. 2000, TERM; 
ACTIVITY, , XX (1) .EQ.500; 

ACTIVITY, ,XX(1)  .GT.500,ZAAL; 
ACTIVITY, , XX (1) . LT. 500, TWOS; 

CLK5  ASSIGN,XX{2)=TNOW,l; 

ACTIVITY; 

ZAAL  GOON,l; 

ACTIVITY, , ,TWOS; 

} 

Q6  QUEUE (6),,,; 

ACTIVITY(l) /6,ATRIB{6) ; 

ASN6  ASSIGN, XX{1)=XX(1)+1,1; 

ACTIVITY, , XX (1) .GE.2000,TERM; 
ACTIVITY, ,XX(1) .EQ.500; 
ACTIVITY,,XX(1)  .GT.500,ZAAM; 
ACTIVITY, ,XX{1) .LT. 500, TWOS; 

CLK6  ASSIGN,XX(2)=TNOW, 1; 

ACTIVITY; 

ZAAM  GOON,l; 

ACTIVITY, , , TWOS ; 

/ 

ASN7  ASSIGN, ATRIB(8)=TNOW-ATRIB(l) ,1; 
ACTIVITY; 

COLl  COLCT(l) ,ATRIB{8) , RESPONSE  TIME, ,1; 
ACTIVITY, ,  ,ASNA; 

f 

TERM  TERMINATE,!; 

END; 

FIN; 


Service  center  5 

Count  event 
End  replication 
Start  stats 
Collect  stats 
Warm-up  period 
Record  stat  start  time 


Service  center  6 

Count  event 
End  replication 
Start  stats 
Collect  stats 
Warm-up  period 
Record  stat  start  time 


Sojourn  time 
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A.8  FORTRAN  Subroutine  for  Model  O2 


•k'k'k'k'k'ii'k-k'k'k-k'k'k'k'k’k-k-k-k-k-k-k-^-k-it-ifk-k-k-k-k-k-k-if-kic-k-k-k-kifk-k-k-k'k'k'k-k-k-ii-k 

*  Main  program  * 

PROGRAM  MAIN 
DIMENSION  NSET(5000) 

INCLUDE  'PARAM.INC 

COMMON/SCOMl/ATRIB(MATRB) ,  DD(MEQT),  DDL(MEQT),  DTNOW,  II,  MFA, 
1MST0P,NCLNR,  NCRDR,  NPRNT,  NNRUN,  NNSET,  NTAPE,  SS (MEQT) , 

2SSL (MEQT) ,TNEXT,  TNOW,  XX (MMXXV) 

COMMON  QSET(5000) 

EQUIVALENCE  (NSET (1) , QSET (1) ) 

NNSET=5000 

NCRDR=5 

NPRNT=6 

NTAPE=7 

OPEN (UNIT=NCRDR, FILE= ’ fort . 5 ' ) 

OPEN (UNIT=NPRNT, FILE= ' fort . 6 ’ ) 

CALL  SLAM 
STOP 

END 

*  Subroutine  to  calculate  stats  and  controls  at  end  of  each  run  * 

*  Subroutine  variable  definitions 

*  Y{NNRUN)  =  average  return  time  for  each  replication  (response) 

*  RATE (NNRUN)  =  average  event  rate 

*  WAIT (NNRUN, J)  =  average  wait  at  J'th  station  for  each  replication 

*  Number  =  number  of  replications 

SUBROUTINE  OTPUT 
INCLUDE  'PARAM.INC 

COMMON/SCOM1/ATRIB(MATRB) ,  DD(MEQT),  DDL (MEQT),  DTNOW,  II,  MFA, 
1MST0P,NCLNR,  NCRDR,  NPRNT,  NNRUN,  NNSET,  NTAPE,  SS(MEQT), 

2SSL (MEQT) ,TNEXT,  TNOW,  XX (MMXXV) 

DOUBLE  PRECISION  Y(IOOO),  WAIT  ( 1000, 5 ) ,  UTILIZ(IOOO) 

INTEGER  matOpen,  mxCreateFull,  matClose 
INTEGER  mxGetPr,  MatPut  Matrix 
INTEGER  b,  d,  FP,  STAT 
INTEGER  J,  NUMBER 

NUMBER  =  1000 ’ 

*  Calculate  system  response  characteristics  * 

-k-k-ir-k-k-ir'k'k'k-k'k-k-k-k-k'ifk-k-k-k-k-ififk-ie-k'k-k-k-k'k'k-k'k-k^-k'k'k-k'k'k-k'k-k-k-k-k-if-k-k-k’k-k-k’k'k'k'k 


*  Calculate  average  response  time  for  each  replication 
Y (NNRUN)  =  CCAVG(l) 

UTILIZ (NNRUN)  =  (CCAVG (14) *CCNUM(14) ) / (TNOW-XX (2) ) 
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*  Output  response  data  in  Matlab  format 

IF  (NNRUN  ..GE.  NUMBER)  THEN 

FP  =  matOpen (' external .mat ' ,  'w') 
b=  mxCreateFull {NUMBER, 1,0) 

call  mxCopyReal8ToPtr (Y,  mxGetPr(b),  NUMBER) 
call  mxSetName(b,  ' extcontrol ’ ) 
stat  =  matPutMatrix (fp,  b) 

call  mxCopyReal8ToPtr (UTILIZ,  mxGetPr(b),  NUMBER) 
call  mxSetName(b,  'extutiliz') 
stat  =  matPutMatrix (fp,  b) 
call  mxFreeMatrix (b) 


d=  mxCreateFull (NUMBER, 5, 0) 

call  mxCopyReal8ToPtr (WAIT,  mxGetPr(d),  5*NUMBER) 

call  mxSetName(d,  'extwait') 

stat  =  matPutMatrix (fp,  d) 

call  mxFreeMatrix (d) 

stat  =  matClose(fp) 

END  IF 


RETURN 

END 
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Appendix  B:  MA  TLAB  Script  Files 
B.  1  MATLAB  Main  Program  for  Analysis  of  Simulation  Data 


%  Computes  and  compares  controlled  responses  from  SLAM  data 

% 

%  Load  data 

% 

load  model; 
load  external; 
load  analytical . out  -ascii; 
load  cntl; 
load  mu; 
load  extmu; 

% 

%  Set  response  and  analytical 
%  analyzed 

% 

response= [sojourn  utiliz] ; 
ext= [extcontrol  extutiliz] ; 
int=[stdwork  stdroute] ; 

[runs, number] =size (response)  ; 
for  q=l: number 

resp=response ( : , q) ; 
analytcv=analytical ( : ,  q) ; 
extcontrol cv=ext ( : ,  q)  ; 

% 

%  Calculate  results  for  10  reps/100  exps,  then  20  reps/50  exps 
% 

for  p=l:2 

reps=10*p; 

exps=100/p; 

% 

%  Perform  analysis  on  {reps}  replications  with  {exps}  experiments 
% 

for  lc=l:exps 

y=resp { (k-1) *reps+l ; k*reps, : ) ; 
intcv  =int { ( k-1 ) *reps+l : k*reps, : ) ; 
analyt=analytcv( (k-1) *reps+l ; k*reps, : ) ; 
extcv=extcontrolcv( (k-1) *reps+l : k*reps, : ) ; 

% 

%  Compute  uncontrolled  response 

% 

uncntl (k, : ) =uncntlst (y,mu (q) ) ; 

% 

%  Compute  controlled  responses 

% 

anltcntl (k, : ) =control (y,mu (q) , analyt, extmu (q) ) ; 
ext cntl (k, : ) =control (y,mu (q) , extcv, extmu (q) ) ; 
int_temp=control (y,mu (q) , intcv, 0.0) ; 
total=size (int) ; 
for  r=l:  {2"'total  (1,2)  )-l 
row= (r-1) *exps+k; 
intcntl ( row, : ) =int_temp ( r, : ) ; 
end  Sfor  loop-r 

end  %for  loop-k  (all  experiments  of  10/20  reps) 

% 


%  Internal  cv's  (stdwork,  stdroute) 

%  External  cv's  (extcontrol  extutiliz) 
%  Analytical  cv' s 

%  t-statistics  and  binary  matrices 
%  expected  values  of  response  stats 
%  expected  values  of  analytical  model 

and  external  control  variates  to  be 
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%  Calculate  average  results  for  10/20  reps 

% 

%  Calculate  average  responses 

% 

if  p==l 

resultlO (q, : ) =mean (uncntl) ; 

% 

%  Calculate  controlled  responses 

% 

redlO_int=result (resultlO (q,  : ) , q,p, intcntl, total (1,2) ) ; 
redlO_anal=result (resultlO (q, : ) , q,p, anltcntl, 1) ; 
redlO_ext=result (resultlO (q,  : ) ,q,p, extcntl, 1) ; 
if  q==l 

redl0= [ redl0_int ; redlO_anal ; redl0_ext ] ; 
else 

redl 0= [ redl 0 ; redl 0_int ; redl 0_anal ; redl 0_ext ] ; 
end  %if-q 

clear  uncntl  anltcntl  extcntl  intcntl  int_temp 

% 

else 

% 

%  Calculate  averages  for  20  reps 
% 

result20 (q, : ) =mean (uncntl) ; 

% 

%  Calculate  controlled  responses 
% 

red20_int=result ( result20 (q,  : ) ,  q,p, intcntl, total (1,2) ) ; 
red20_anal=result (result20 (q,  : ) , q, p, anltcntl, 1) ; 
red20_ext=result (result20 (q,  : ) , q,p, extcntl, 1) ; 
if  q==l 

red2  0= [ r ed2  0_int ; red2  0_anal ; red2  0_ext ] ; 
else 

red2  0= [ red2  0 ; red2  0_int ; red2  0_anal ; red2  0_ext ] ; 
end  %if-q 

end  %if-p  (100  or  50  experiments) 
end  %for  loop-p  (10  or  20  runs) 
end  %for  loop-q  (response  under  study) 

% 

%  Save  results 

% 

save  resultlO. out  resultlO  -ascii 
save  redlO.out  redlO  -ascii 
save  result20.out  result20  -ascii 
save  red20.out  red20  -ascii 
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B.2  MATLAB  Function  for  Calculating  Controlled  Responses 


function  z  =  control (y, mu, x, extmu) 
load  cntl; 

% 

%  Returns  controlled  responses  and  associated  statistics 
%  Input : 

%  y  =  uncontrolled  response 

%  mu  =  known  mean  of  y 

%  X  =  control  variates 

%  extmu  =  known  mean(s)  of  x 

%  Output : 

%  z  =  [controlled  mean,  variance,  confidence  interval  width, 

%  coverage  (0  or  1),  mean  square  error  from  mu] 

% 

[len,wid]=size (x) ; 

% 

%  Load  appropriate  binary  matrix  to  compute 
%  all  possible  combinations  of  cv's 

% 

bin=binary (wid) ; 

% 

%  compute  z  for  all  possible  cv  combinations 
% 

comb= (2^wid) -1; 
e=ones (len, 1) ; 
x=x- extmu; 
for  j=l:comb 

X=[e  x( : ,bin( j, : ) ) ] ; 
beta=X\y; 
z(j,l)=beta(l,l) ; 

[l,w] =si2e (X) ; 
sigma=inv (X ' *X) ; 

SSE=y ’ *y-beta ’ *X ' *y; 

MSE=(1/ (l-w) ) *SSE; 

z ( j , 2 ) =sigma (1,1) *MSE; 

z  { j, 3) =2*t (len-w, : ) *sqrt (z ( j,2) ) ; 

ul=2{j,l)+.5*z(j,3)  ; 

ll=z(j,l)-.5*z(j,3)  ; 

z ( j ,  4 ) =mu<=ul  &  mu>=ll; 

z(j,5)  =  (z(j,l)-mu)''2; 

end 
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dp  dp  df> 


B.3  MATLAB  Function  for  Computing  Uncontrolled  Response  and  Statistics 


function  z=uncntlst (y,mu) 
load  cntl; 

% 

%  Computes  uncontrolled  response  and  statistics 
%  Input: 

%  y  =  uncontrolled  response 

%  mu  =  true  mean 

%  Output : 

z  =  [mean,  variance,  confidence  interval  width,  coverage 
(0  or  1),  mean  square  error  from  true  mean] 


z  (1)  =mean  (y)  ; 

z  (2)  =  (std (y) ^2) /length (y)  ; 

z (3) =2*t (length (y) -1, : ) *sqrt (z (2) ) ; 

u1=2(1)+.5*z{3) ; 

ll=z(l)-.5*z(3); 

z ( 4 ) =mu<=ul  &  mu>=ll ; 

z(5)  =  (z(l)-mu)'^2; 
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B.  4  MA  TLAB  Function  to  Average  Statistics  over  all  Experiments 

function  z=result (y,q,n,x,w) 

% 

%  Calculate  average  controlled  responses  for  100/n  experiments 

% 

%  Input : 

%  y  =  uncontrolled  responses 

%  q  =  response  statistic  indicator 

%  n  =  counter  to  determine  nximber  of  experiments 

%  X  =  controlled  responses 

%  w  =  number  of  control  variates 
%  Output : 

%  z  =  [q,  variance,  confidence  interval  width,  coverage,  MSE, 

%  variance  reduction  percent,  confidence  interval  percent] 

% 

exps=100/n; 
for  m=l:  (2-'w)-l 
z(m, l)=q; 

z (m, 2 ) =mean (x ( (m-1 ) *exps+l :m*exps, 2 ) ) ; 
z (m,  4) =mean (x ( (m-1) *exps+l :m*exps, 3) ) ; 
z (m,  6:7) =mean (x ( (m-1 ) *exps+l :m*exps,  4:5)); 
z(m,3)=(l-(z(m,2)/y(l,2) ))*100; 
z(m, 5)=(l-(z(m, 4)/y(l,3) ) )*100; 
z  (m,  8) =mean (x ( (m-1 ) *exps+l :m*exps, 1) ) ; 

end 
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B.5  MA  TLAB  Function  for  Returning  Appropriate  Binary  Matrix 

function  z=binary (wid) 

% 

%  Returns  appropriate  "binary"  file  for  number  of  control  variates 
%  (For  example  if  number  of  cv' s  is  3,  then  appropriate  binary 
%  matrix  is:  [001;  010;  100;  011;  101;  110;  111]  ) 

% 

%  Input : 

%  wid  =  number  of  control  variates 
%  Output : 

%  z  =  [binary  matrix] 


if  wid==l 
bin=binl; 
el seif  wid==2 
bin=bin2 ; 
elseif  wid==3 
bin=bin3 ; 
elseif  wid==4 
bin=bin4 ; 
elseif  wid==5 
bin=bin5; 
elseif  wid==6 
bin=bin6; 
if  wid==7 
bin=bin7 ; 
elseif  wid==8 
bin=bin8; 
elseif  wid==9 
bin=bin9; 
elseif  wid==10 
bin=binl0; 
elseif  wid==ll 
bin=binll; 
end 
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